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A  B  ST RACT 


This  Semiannual  Technical  Summary  describes  the  Lincoln  Laboratory 
Vela  Uniform  programmer-  the  period -4- October -lfraO-ihrough  it  March- 
t?SiT^  Section  I  describes  progress  in  the  development  of  a  prototype 
Seismic  Data  Center.  This  development  is  on  schedule,  and  a  first 
version  of  the  prototype  will  be  operational  during  the  next  six  months. 
Section  II  describes  a  series  of  activities  in  the  development  and  imple¬ 
mentation  of  the  seismic  processing  techniques  that  will  be  utilized  in 
the  prototype  system.  Section  III  describes  a  series  of  investigations 
in  General  Seismological  Research. 
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. SUMMARY 


This  la  the  thirty-fourth  Semiannual  Technical  Summary  (SATS)  report  describing  the 
activities  of  Lincoln  Laboratory  funded  under  Project  Vela  Uniform.  This  report  cove  ra  the 
period  1  October  1980  through  31  March  1981.  Project  Vela  Uniform  Is  a  program  of  research 
Into  the  discrimination  between  earthqt*  kes  and  nuclear  explosions  by  seismic  means.  An 
Important  recent  emphasis  of  the  project  la  in  the  development  of  data-handling  and  analysis 
techniques  that  may  be  appropriate  for  the  monitoring  of  a  potential  Comprehensive  Test  Ban 
Treaty.  The  FY  1981  program  consists  of  the  development  of  a  prototype  Seismic  Data  Center 
(SDC),  which  will  fulfill  U.  S.  obligations  that  may  be  Incurred  under  such  a  Treaty,  and  which 
will  also  provide  a  data  resource  to  the  selsmological  research  community.  The  program  also 
contains  an  element  of  seismic  research,  with  emphasis  on  those  areas  directly  related  to  the 
operations  of  the  SDC. 

Substantial  progress  has  been  made  in  the  development  of  the  first  version  of  the  prototype 
SDC.  In  the  absence  of  a  Local  Computer  Network,  which  will  be  delivered  during  the  next 
quarter,  the  Seismic  Analysis  Station  (SAS)  and  database  subsystems  are  being  developed  as 
stand-alone  systems.  Seismic  processing  techniques  are  being  developed  on  our  general- 
purpose  computers.  We  have  now  demonstrated  the  capability  to  read  Input  seismic  data, 
carry  out  signal  detection  and  the  automatic  extraction  of  a  variety  of  signal  parameters,  and 
display  the  detected  waveforms  for  analyst  review.  A  variety  of  additional  analyst  software 
tools  will  be  installed  soon  on  the  SAS.  The  waveform  and  parameter  database  systems  are 
being  developed  on  the  SAS  as  local  facilities,  and  they  will  be  expanded  to  system-wide  data¬ 
bases  at  a  later  stage.  A  Remote  Seismic  Terminal,  capable  of  limited  waveform  display 
and  analysis  and  remote  interaction  with  the  SDC  over  dial-up  lines,  has  been  demonstrated. 

An  important  effort  during  this  period  has  been  the  development  and  implementation  of  a 
variety  of  seismic  processing  algorithms  for  automatic  signal  analysis,  and  analyst  interactive 
review,  A  package  of  programs  for  automatic  analysis  of  incoming  data  has  been  installed  on 
our  VAX  li/780  computer.  These  include  programs  for  signal  detection,  onset  time  measure¬ 
ment.  amplitude  and  period  measurement,  and  discrimination  between  local  and  teleseismic 
events.  In  the  prototype  SDC,  provision  is  made  for  analyst  review  and  modification  of  each 
of  these  measurements.  The  output  sets  of  arrival  parameters  are  then  merged  with  any  ex¬ 
ternal  arrival  data  into  a  local  parameter  database.  This  is  follow ed  by  automatic  association 
into  a  preliminary  event  list,  analyst  review  of  the  output  using  available  waveform  data,  and 
final  event  bulletin  preparation.  We  are  continuing  to  evaluate  and  develop  a  series  of  addi¬ 
tional  algorithms  for  phase  Identification,  azimuth  measurement,  and  association  of  long- 
period  arrivals  with  events  determined  from  short-period  data. 

Other  seismological  research  has  continued  in  several  areas.  Previous  work  on  anom¬ 
alous  surface-wave  excitation  by  presumed  explosions  in  the  Eastern  Kazakh  region  has  been 
extended.  A  variety  of  evidence  is  given  in  support  of  the  hypothesis  that  thrust  faulting  ac¬ 
companies  these  unusual  events,  and  that  a;  least  some  observed  tectonic  strain  release  by 
Nevada  Test  Site  events  can  also  be  attributed  to  this  type  of  faulting.  Other  studies  include: 
an  inversion  of  Rayleigh-wave  phase  and  group  mantle  structure  in  the  Middle  East  and  South 
Central  Asia;  a  method  for  the  evaluation  of  mode  eigenfunctions  for  P  and  SV  motion  in  a 
layered  medium,  which  is  numerically  stable;  and  an  application  of  this  technique  to  the  study 
of  the  excitation  of  the  Lg  phase  as  a  function  of  source  depth.  A  generalization  of  the  method 
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of  centroid  locatioh  to  the  estimation  of  sourer  parameters  (particularly  moment  and  hypoc enter) 
for  events  close  in  time  is  described  and  applied  to  the  resolution  of  explosion  and  earthquake 
sources,  and  to  the  modeling  of  finite  seismic  sources.  In  another  investigation,  the  effect  of 
a  boundary  on  the  radiation  from  seismic  sources  is  evaluated,  with  particular  emphasis  on 
contributions  to  the  moment  tensor.  It  is  shown  that  assumptions  about  the  boundary  are  nec¬ 
essary  if  a  proper  interpretation  of  the  moment  tensor  Is  to  be  obtained.  The  application 
of  finite  difference  methods  to  the  calculation  of  normal-mode  eigenfrequencles  and  eigen¬ 
functions  so  that  seismogram  synthesis  may  be  performed  at  higher  frequencies  is  discussed. 

M.  A.  Chinnery 
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SEISMIC  DISCRIMINATION 

t.  SEISMIC  DATA  CENTER  DEVELOPMENT 

A.  PROGRESS  IN  SEISMIC  DATA  CENTER  (SDC)  DEVELOPMENT 

During  this  report  period,  development  has  continued  on  the  various  subsystems  of  the 
1 

prototype  SDC.  Particular  emphasis  is  being  placed  on  the  development  of  the  SeiBmic  Analysis 
Station  (SAS),  the  waveform  database,  and  the  parameter  database.  Hardware  acquisition  for 
these  subsystems  is  now  complete.  The  Local  Computer  NetWork,  which  will  link  these  sub¬ 
systems,  is  being  built  by  a  subcontractor  and  is  scheduled  for  delivery  within  a  few  months. 

The  first  operational  version  of  the  SDC  is  being  built  up  around  a  combination  of  the  SAS 
and  a  general-purpose  (VAX  11/780)  computer,  with  local  waveform  and  parameter  databases. 
This  early  version  will  provide  both  a  test  bed  for  the  refinement  of  the  man-machine  interface 
and  will  also  help  to  define  the  structure  of  the  databases  for  implementation  on  a  system-wide 
basis  early  in  FY  82.  The  VAX  also  performs  all  seismic  processing  tasks,  and  is  therefore 
acting  as  a  test  bed  for  the  Automatic  Processing  Subsystem,1 

Initial  software  implementation  on  the  VAX  and  SAS  is  essentially  complete.  Seismic  data 
received  as  part  of  the  International  Data  Collection  Experiment  (IDCE)  can  now  be  processed 

using  a  variety  of  algorithms  (see  Sec.  II  of  this  report).  The  local  waveform  database  is  a  de- 

2 

velopment  of  the  "gram  file"  concept,  and  the  local  database  is  an  improvement  on  the  "Mod  I 

system."3  Each  of  these  will  be  available  to  the  analyst  using  the  SAS.  Display  software  for 

2 

the  SAS  has  been  based  on  a  "display"  system  evolved  previously.  Each  of  these  components 
is  described  in  the  contributions  that  follow. 

During  the  past  six  months,  substantial  effort  has  also  been  expended  on  the  development 
of  a  Remote  Seismic  Terminal  (RST),  which  will  form  a  low-cost  remote  facility  for  limited 
waveform  analysis  and  interaction  with  the  SDC.  This  unit  was  demonstrated  to  the  UN  Com¬ 
mittee  on  Disarmament  in  February  1961. 

Specifications  for  the  Communications  Interface  Subsystem,  which  will  link  the  SDC  to  in¬ 
coming  data  streams,  are  being  formulated  on  the  basis  of  discussions  with  DOE.  This  will  be 
built  by  a  subcontractor  during  the  next  six  months. 

We  continue  to  interact  with  DOE,  DARPA,  and  various  DARPA  subcontractors  to  ensure 
that  our  approach  to  the  prototype  SDC  is  state-of-the-art, 

M.  A.  Chinnery 
A.  G.  Gann 

B.  SEISMIC  PROCESSING  IN  THE  PROTOTYPE  SDC 

This  section  describes  the  seismic  processing  steps  that  will  be  performed  in  the  Version  1 
of  the  Seismic  Data  Center  prototype,  SDC-V1,  The  goal  of  SDC-V1  is  to  provide  many  of  the 
capabilities  of  the  SDC  as  quickly  as  possible  so  that 'operational  experience  with  SDC-V1  can 
lead  to  rapid  improvements  in  later  versions.  The  capabilities  of  SDC-V1  will  include: 

(1)  Automatic  detection  of  arrivals  using  waveform  data  from  various 
sources. 

(2)  Analyst  review  and  refinement  of  automatic  detections  using  interactive 
graphics. 
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(3)  Combining  automatic  detections  and  external  arrival  times  to  produce 
an  arrival  database, 

(4)  Combining  external  event  bulletin  Information  Into  a  parameter  database, 

(5)  Automatic  association  of  arrivals  to  produce  a  preliminary  bulletin, 

(6)  Analyst  review  and  refinement  c  f  un  preliminary  bulletin. 

(?)  Production  of  final-event  bulletins  based  on  several  event-selection 
criteria. 

Routine  analysis  will  have  the  following  major  steps; 


DP 

Detection  Processing 

ABAR 

Arrival  Based  Analyst  Review 

MERGE 

Merge  Parametric  Data  Sources 

AA 

Automatic  Association 

EBAR 

Event  Based  Analyst  Review 

BUL 

Bulletin  Preparation 

Normally,  processing  will  proceed  from  one  stage  to  the  next,  but  circumstances  occasion¬ 
ally  will  require  that  one  or  more  stages  must  be  redone. 

The  following  subsections  briefly  describe  the  processing  and  software  involved  in  each 
stage. 


1.  DP  —  Detection  Processing 

Inputs;  Continuous  waveform  data  stored  in  GRAMFILE  format.  Waveform  data  are  cur¬ 
rently  received  on  magnetic  tape  in  various  formats.  Programs  which  read  the  tapes  and  store 
the  waveforms  in  a  GRAMFILE  database  are  available. 

Outputs;  A  GRAMFILE  database  which  contains  waveform  segments  of  the  interval  around 
the  arrival  time,  and  a  time-ordered  list  of  detected  arrival  times  and  related  parameters. 

Processing;  The  detection  processing  is  done  in  two  steps.  First,  a  detection  algorithm 
is  run  on  the  continuous  waveforms  to  detect  arrivals  and  produce  waveform  segments.  A 
second  process  is  run  on  the  segmented  data  to  refire  the  onset  time  and  measure  other 
parameters. 

2.  ABAR  —  Arrival  Based  Analyst  Review 

Inputs :  The  outputs  of  the  previous  step. 

Outputs ;  The  arrival  list  and  waveform  database  are  revised  to  reflect  improvements  made 
by  the  analyst. 

Processing;  Using  the  DISPLAY  command,  the  analyst  inspects  the  waveform  data  corre¬ 
sponding  to  each  pick,  one  at  a  time,  and  makes  any  adjustments  and  additional  measurements 
that  might  be  required.  Waveform  data  which  have  been  provided  to  the  SDC  in  segmented  form, 
such  as  from  SRO  stations,  can  be  inspected  in  this  way  also.  Once  the  analyst  is  finished,  the 
second  stage  of  the  DP  is  run  to  redo  any  automatic  measurements  as  necessary. 

3.  MERGE  —  Merge  Parametric  Data  Sources 

Inputs;  Internal  arrivals  from  the  previous  step  and  arrival  data  from  external  sources 
such  as  the  WMO/GTS, 
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Outputs :  A  uniform  parametric  database  for  a  given  time  interval.  Thia  database  will  con* 
tain  at  least  a  time-ordered  arrival  list,  and  a  time-ordered  event  list  if  event  information  is 
available  from  external  sources.  Currently,  this  database  is  in  MODI  PARAMETERFILE 
format. 

Processing:  To  date,  the  parametric  data  created  by  the  previous  steps  have  been  kept  in 
a  form  convenient  for  the  detection  and  arrival  review  process,  in  separate  databases  for  each 
station  for  example.  These  data  will  be  merged  with  any  external  arrival  or  event  information 
to  produce  a  parametric  database.  External  data  may  require  extensive  preprocessing  and 
editing.  The  sorting  and  editing  tasks  are  performed  using  standard  UNIX  commands,  and 
EDPF  -  the  PARAMETERFILE  editor. 

4.  AA  -  Automatic  Association 

Inputs;  The  parametric  database  produced  in  the  previous  step. 

Outputs:  A  parametric  database  augmented  with  event  locations,  and  GRAMFILE  databases 
for  each  event  for  use  with  DISPLAY  during  EBAR. 

Processing:  An  Automatic  Association  (AA)  algorithm  is  run  on  the  parametric  database. 

AA  searches  the  arrival  list  and  identifies  group  arrivals  that  belong  to  an  esrthquake.  The 
arrivals  are  then  used  to  locate  the  event,  and  the  event  is  added  to  the  database.  This  process 
is  repeated  for  the  entire  arrival  list.  The  resulting  database  can  be  used  to  produce  a  prelim¬ 
inary  event  bulletin. 

After  association,  the  waveform  data  are  grouped  by  event  for  uBe  in  the  following  step. 

5.  EBAR  —  Event  Based  Analyst  Review 

Inputs:  Outputs  from  the  previous  step. 

Outputs:  A  refined  parameter  database. 

Processing:  An  analyst  is  again  given  an  opportunity  to  review  the  results  of  automatic 
processing,  using  DISPLAY  in  conjunction  with  the  parameter  database  editor,  DEV.  However, 
at  this  point  he  has  much  more  information  to  go  on.  namely  the  event  hypotheses  produced  by 
AA.  The  analyst  may  need  to  do  many  operations  on  the  waveform  and  parametric  data.  For 
example,  to  refine  an  event  location  an  analyst  may  need  to  display  the  waveforms  corresponding 
to  a  subset  of  the  arrivals  associated  with  this  event,  adjust  picks  or  other  parameters,  alter 
associations,  and  relocate  the  event.  This  process  may  have  to  be  Iterated  several  times  for 
each  event. 

6.  BUL  —  Bulletin  Preparation 

Inputs:  The  parametric  database  from  the  previous  step. 

Outputs:  Various  event  bulletins, 

Processing;  Once  the  parameter  database  has  been  processed  to  the  satisfaction  of  the 
analyst,  event  bulletins  can  be  produced  as  needed.  This  task  will  be  done  by  a  senior  analyst 
who  is  responsible  for  bulletin  preparation  and  ensuring  that  the  parameter  database  is  complete, 
consistent,  and  that  its  analysis  is  complete. 


K.  R.  Anderson 


C.  WAVEFORM  DISPLAY  ON  THE  SAS 

The  SAS  is  the  principal  user  terminal  for  performing  the  operational  and  research  functions 
of  the  SDC.  The  current  prototype  consists  of  a  PDP-ii/44  computer,  Megatek  7000  graphics 
display,  two  Datagraphix  132 -A  alphanumeric  displays,  disks,  and  tapes.  This  subsystem  will 
be  linked  to  the  rest  of  the  SDC  subsystem  through  the  local  computer  network  when  it  is  in- 
stalled.  The  current  prototype  is  being  used  to  evaluate  display  techniques  for  waveform  data, 
the  integration  of  waveform  and  parameter  data  processing,  and  the  development  of  enhanced 
software  for  the  SAS  subsystem  itself, 

A  program  called  DISPLAY  was  developed  several  years  ago  on  the  PDP-li/50  to  allow 
seismologists  to  Interactively  analyse  a  waveform  database  on  the  Tektronix  4014  graphics 
terminal.  DISPLAY  plots  a  subset  of  the  waveforms  in  a  database,  and  then  allows  the  user  to 
interactively  set  and  delete  markers,  change  vertical  and  horizontal  scaling,  post  amplitude 
and  time  scales,  and  move  the  waveforms  left  and  right.  We  have  converted  DISPLAY  to  run 
on  the  prototype  SAS.  The  graphics  primitives  originally  in  DISPLAY  were  replaced  by  calls 
to  Graphpac,  a  faster  set  of  graphics  routines  which  originally  were  also  used  on  the  Tektronix. 
Conversion  of  DISPLAY  to  the  Megatek  was  completed  by  using  a  set  of  device  specific  graphics 
routines  developed  by  Purdue  University,  As  it  now  stands,  DISPLAY  calls  Graphpac  which, 
in  turn,  calls  the  Purdue  routines  to  display  information  on  the  Megatek  graphics  terminal.  All 
the  options  available  in  the  original  DISPLAY  program  designed  for  the  Tektronix  4014  graphics 
terminal  can  now  be  used  with  the  Megatek  7000  graphics  terminal  en  the  PDP-11/44. 

Complementing  this  approach,  we  are  working  on  a  variety  of  performance  enhancements. 
Here,  software  at  the  device  specific,  function  primitive,  and  system  levels  will  be  designed 
and  implemented  to  facilitate  new  waveform  manipulations,  for  which  the  original  software  is 
too  deficient  structural'y  to  be  modified  appropriately. 

Intrinsic  in  this  dual  approach  .s  a  layered,  or  fine-grained  approach  to  software  develop¬ 
ment.  Conceptually,  three  levels  or  layers  of  software  are  involved;  (1)  a  command  level; 

(2)  the  function  primitive  level;  and  (3)  the  device  specific  or  system  levels.  The  command 
level  is  a  user  interface  accepting  input  from  the  user  and  dividing  his  commands  into  more 
elementary  function  primitives  related  to  display  and  database  tasks.  The  function  primitives 
retrieve  and  process  data  necessary  for  display  tasks  and  create  various  device  specific  and 
system  level  function  calls  necessary  to  drive  the  graphics  and  other  hardware. 

This  layered  structure  underlying  the  software  permits  separation  of  the  user  interface 
from  device  specific  and  lower-level  data-handling  considerations.  As  later  generations  of 
graphics  and  device  specific  primitives  emerge,  permitting,  for  example,  increasingly  com¬ 
plex  waveform  scrolling,  overlaying  and  geometric  transformation  functions,  they  may  be  im¬ 
plemented  in  terms  of  existing  or  slightly  modified  upper  layers  operating  with  a  new  bottom 
level  to  the  system.  As  experience  is  gained  with  command  languages  and  user  interfaces,  in¬ 
creasingly  facile  ones  can  be  implemented  while  retaining  existing  lower  levels  to  the  system. 

At  any  time,  as  at  present,  a  working  system  is  available  which  can  be  utilized  for  seismic 
analysis  and  as  a  laboratory  for  discovering  improved  graphics  techniques  for  seismic  analysis. 

P.  T.  C rames 
J.  Hirsch 


D.  WAVEFORM  DATABASE  SUBSYSTEM 

The  Database  Subsystem  of  the  SDC  consists  of  two  conceptually  separate  parts,  the  wave¬ 
form  database  and  the  parameter  database.  The  waveform  database  must  receive  the  telemetered 
waveform  data  from  the  communications  interface  subsystem,  store  it  in  the  on-line  working 
set.  archive  it  in  permanent  storage,  and  process  retrieval  requests  for  specific  blocks  of 
waveform  data.  The  parameter  database  receives  parametric  data  from  various  subsystems 
in  the  SDC,  updates  the  database,  and  processes  retrieval  requests  for  parametric  data  for 
other  subsystems  in  the  SDC. 

This  section  describes  a  preliminary  specification  for  the  waveform  database  component 
of  the  SDC.  Two  systems  are  mentioned:  the  initial  prototype,  which  should  be  completed  by 
October  1981;  and  the  final  prototype,  which  is  projected  for  FY  1982.  The  initial  prototype  is 
a  simple  system  that  stores  only  waveform  data;  parametric  data  will  be  stored  locally  at  the 
Seismic  Analyst  Station.  The  final  prototype  will  integrate  both  waveform  and  parametric  data 
storage  and  retrieval. 


1.  Waveform  Database  Sizing  Requirements 

The  final  prototype  must  be  able  to  absorb  incoming  data  at  a  rate  of  200  kbps.  This  figure 
is  double  the  expected  100-kbps  nominal  data  rate  from  30  National  Seismic  Station  (NSS)  in¬ 
struments.  Doubling  the  expected  rate  allows  for  retransmission  after  perious  of  Data  Center 
shutdown  or  communications  failure.  Also,  it  will  allow  moderate  increases  in  the  nominal 
data  rates  caused  by  additional  NSSs  and  other  data  sources,  such  as  Seismic  Research  Obser¬ 
vatory  (SRO)  data  on  magnetic  tape.  The  design  should  permit  expansion  beyond  this  data  rate 
since  additional  NSSs  and  other  data  sources  may  be  added  to  the  system.  However,  the  desire 
for  expandability  should  not  drive  the  design  to  unrealistically  ambitious  aims. 

The  requirements  for  the  initial  prototype  will  be  much  less  demanding.  The  data  source 
will  be  five  NSSs  of  the  Regional  Seismic  Test  Network  (RSTN).  Each  station  provides  roughly 
2.4  kbps,  so  the  incoming  data  rate  will  be  about  12  kbps.  There  are  no  initial  reliability  re¬ 
quirements,  so  a  guideline  dat-  rate  of  20  kbps  should  allow  for  adequate  expansion. 

Both  the  initial  and  final  prototypes  should  be  able  to  retain  all  incoming  data  on-line  during 
a  5-day  bulletin  preparation  period.  The  resulting  capacity  requirements  for  urorocessed  data 
during  various  time  periods  are  displayed  in  Table  1-1.  Note  that  tl  se  calculations  assume  no 
data  conversion  or  data  compression  since  disk  formats  have  not  yet  been  defined.  In  addition, 
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TABLE  1-1 

WAVEFORM  SIZING  REQUIREMENTS 


Data  Rates 

Initial  Prototype 

Final  Prototype 

Kilobits/second 

20.0 

200 

Kilobytes/second 

2.5 

25 

Megabytes/day 

216.0 

2,160 

Megabytes/5  days 

1,080.0 

10,800 

i.d.-  linnilij  ViiLiliiiHilil  rrlililhK ,J| 


space  will  be  required  for  indices  and  other  system -controlled  areas,  and  for  parametric  data 
in  the  final  prototype, 

2.  Reliability  Requirements 

The  following  reliability  requirements  apply  to  the  final  prototype  system.  They  are  not 
required  for  the  initial  prototype. 

In  general,  data  should  continue  to  be  stored  somewhere  within  the  Data  Center  during  rou¬ 
tine  maintenance  and  small-scale  component  failures.  Retransmission  from  external  sources 
should  be  required  only  for  recovery  from  extensive  failures.  This  means  that  either  the  wave¬ 
form  database  must  use  a  non-stop  implementation  or  the  Communications  Interface  must  be 
able  to  retransmit  data.  Details  of  these  approaches  will  be  examined  and  the  choice  of  approach 
will  be  made  during  the  design  of  the  final  prototype. 

Authentication  of  the  incoming  data  is  performed  by  the  Communications  Interface.  There¬ 
fore,  all  data  reaching  the  waveform  database  are  assumed  to  be  valid.  Some  data  may  be  re¬ 
transmitted  to  the  Data  Center  after  failure  of  the  data  network  supplying  the  data,  or  complete 
failure  of  the  Communications  Interface.  Multiple  transmissions  of  data  should  match  the  orig¬ 
inal  transmission.  Mismatches  should  be  detected  and  some  technique,  perhaps  manual,  should 
be  employed  to  determine  what  data  are  valid. 

3.  The  Initial  Prototype 

The  goal  of  the  initial  waveform  database  prototype  effort  is  to  quicKly  develop  an  imple¬ 
mentation  of  the  waveform  database  that  supplies  the  functionality  essential  to  the  analyst's 
work.  This  approach  will  provide  a  foundation  for  developing  and  refining  functions  that  the 
automatic  processing  tasks  and  analysts  require  the  waveform  database  to  perform.  This  early 
experience  with  an  operational  Data  Center  is  essential  to  determine  what  features  are  really 
required,  and  will  guide  the  design  of  the  final  prototype. 

The  initial  prototype  waveform  database  will  be  implemented  using  the  UNIX^  operating 
system  on  a  single  DEC  VAX  li/780  processor  with  4  Mbytes  of  primary  memory.  The  con¬ 
figuration  will  include  four  600-Mbyte  disk  drives  for  waveform  storage  during  bulletin  prep¬ 
aration.  Received  data  will  be  archived  on  6250-bpi  density  magnetic  tape. 

The  Communications  Interface  will  supply  authenticated  data  to  the  waveform  database  using 
magnetic  tapes.  Alternate  approaches  include  a  connection  to  the  Communications  Interface 
via  a  direct  point-to-point  data  connection  or  the  local  network.  The  final  approach  we  use  will 
depend  on  availability  o'  the  Sytek  Local  Network  Interfaces. 

A  high-speed  point-to-point  data  connection  will  be  used  to  connect  the  Seismic  Analyst 
Station  and  the  waveform  database  since  the  Sytek  Local  Network  Interfaces  will  not  be  available 
soon  enough  to  be  used  for  this  application.  Magnetic  tape  was  ruled  out  for  this  connection 
because  the  extra  manipulation  required  to  transfer  data  via  magnetic  tape  would  limit  the  im¬ 
plementation  to  batch- style  strategies  and  require  undesirable  analyst  action  to  manipulate  the 
tapes. 

Existing  software  wVi  form  the  basis  of  the  waveform  database  implementation.  This  ap¬ 
proach  simplifies  transporting  existing  programs  to  the  VAX,  since  these  programs  already  use 
a  consistent  database  format  and  access  subroutine  library.  It  also  simplifies  implementation 
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of  new  programs  since  our  seismic  research  personnel  are  already  familiar  with  the  existing 
database  formats  and  access  routines.  Also,  the  maturity  of  the  existing  routines  makes  them 
far  more  reliable  than  new  software  would  be. 

The  remaining  development  of  the  initial  waveform  database  involves  enhancements  to  the 
existing  database  software  and  interfacing  the  waveform  database  with  the  Seismic  Analyst 
Station  and  the  Communications  Interface.  The  level  of  effort  required  will  depend  on  how  dif¬ 
ficult  hardware  installation  and  software  performance  and  functional  enhancement  prove  to  be. 

We  are  just  completing  the  hardware  procurement.  The  disk  and  tape  peripherals  have 
arrived  and  are  being  installed.  We  expect  the  intercomputer  links  to  arrive  during  the  xxrst 
week  of  April  1981.  Software  for  these  devices  has  been  developed  at  other  installations  and 
we  are  in  the  process  of  acquiring  copies.  We  expect  installation  to  be  finished  by  the  end  of 
April. 

The  file  system  performance  of  UNIX  should  be  adequate  for  the  2.5-Mbyte/s  design  data 
rates  required  by  the  initial  prototype.  Additional  benchmarks  must  be  done  under  varying 
load  conditions  when  the  large  disks  and  high-density  tapes  have  been  installed  to  determine 
whether  performance  issues  will  require  special  index  structures  not  supplied  by  the  UNIX  file 
system. 


4.  Final  Prototype 

We  are  beginning  the  specification  of  the  final  prototype.  The  design  will  encompass  both 
the  waveform  and  parametric  databases.  Several  critical  design  areas  such  as  user  interface, 
reliability,  extensibility,  and  data  compression  have  been  identified.  We  plan  to  elaborate  and 
refine  the  design  during  the  next  few  months  using  our  experience  with  the  initial  prototype  to 
guide  our  efforts.  The  design  should  be  completed  and  implementation  started  by  October  1981. 

K.  J.  Schroder 

E.  ACCESS  ROUTINES  FOR  THE  PARAMETER  DATABASE 

Work  has  progressed  on  our  parameter  database  access  routines,  a  set  of  file  management 
subroutines  which  enable  the  user  to  create  and  modify  files  of  fixed-length  records  flexibly. 

The  access  routines  have  been  used  by  applications  programs  during  the  Mod  I  project,  and  will 
be  used,  in  enhanced  form,  in  the  development  of  the  SAS.  Briefly,  records  may  contain  pa¬ 
rameters  of  several  different  types;  the  user  describes  the  layout  of  a  typical  record  in  a  C 
structure,  and  then  calls  routines  to  add,  delete,  change,  and  find  records.  To  date,  we  have 
been  able  to  access  x^ese  records  by  two  slow  methods:  physically  sorting  the  file  or  a  desired 
key  and  then  using  binary  searches  to  locate  the  desired  record,  and  sequentially  searching 
through  all  the  records  until  the  desired  key  is  found.  These  access  methods  were  developed 
so  that  we  could  test  the  user  interface  of  the  access  routines  as  soon  as  possible,  without  spe¬ 
cial  regard  for  speed  or  operating  efficiency.  The  sorting  of  records  is  not  efficient  since  we 
can  only  sort  on  one  parameter  at  a  time.  Even  though  fast  searches  can  be  made  on  the  order 
of  log.,N  (where  N  equals  the  number  of  records  in  the  file),  insertions  and  deletions  require 
that,  on  the  average,  N/2  records  be  physically  moved.  The  sequential  searches  are  not  effi¬ 
cient  because,  on  the  average,  N/2  records  must  be  read  before  the  desired  record  is  found. 

We  have  now  developed  a  third  method  of  accessing  records,  B+  trees.  This  method  is 
very  fast,  but  somewhat  more  disk  storage  space  is  required  than  with  the  other  two  methods. 
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B+  trees,  developed  by  R.  Bayer  and  E.  McCreight,  allow  retrieval,  insertion,  and  deletion 
in  time  proportional  to  ing^N,  where  N  is  the  number  of  elements  in  the  index.  In  our  seheme, 
the  value  of  k  depends  on  the  indexed  variable's  type,  as  follows:  k  is.  30  for  variables  typed 
as  double,  60  for  floats,  and  120  for  integers  because  they  are  8-,  4-,  and  2 -byte-long  vari¬ 
ables,  respectively,  in  a  240-byte  fixed-length  record.  Using  floats  as  an  example,  we  can 
access  any  of  60  keys  with  one  disk  access,  3,600  keys  with  two  disk  accesses,  and  216,000  keys 
with  three  disk  accesses.  A  penalty  is  paid  in  terms  of  disk  storage  space,  because  a  B+  tree 
must  be  built  for  each  variable  indexed.  Besides  their  fast  r  andom-access  speed,  B+  trees 
allow  even  faster  sequential  accesses.  Once  the  desired  key  is  located  in  the  tree,  the  next 
higher  or  lower  key  can  be  accessed  directly  because  a  pointer  to  each  is  present. 

P.  T.  Crames 

F.  APPLICATION  OF  THE  MUMPS  SYSTEM  TO  THE  PARAMETER  DATABASE 

A  sample  parametric  database  and  associated  query  capability  were  developed  using  the 
MUMPS  language  and  operating  system  available  on  the  PDP-ll/44,  This  implementation  of 
MUMPS  is  a  stand-alone  system  optimized  for  database  applications.  The  system  includes  the 
operating  system,  data  management  system,  and  query  language.  Sim  e  the  query  language  is 
interpreted,  and  not  compiled,  and  the  basic  primitives  are  aimed  at  string  manipulation,  it 
was  possible  to  develop  a  prototype  system  very  quickly,  in  less  than  two  weeks. 

The  database  was  transferred  from  a  UNIX  system  via  magnetic  tape,  and  loaded  into  a 
MUMPS  database  structure  which  had  been  designed  for  query  access  efficiency  (the  basic  data 
structure  in  MUMPS  is  Multi-way  trees5).  Query  capability  included  selection  of  stations,  time 
intervals,  latitude  and  longitude  ranges,  magnitude,  and  arrival  type.  In  addition,  the  bulletin, 
arrival,  and  event  print  functions  of  the  Mod  I  library  were  duplicated  within  MUMPS. 

Comparison  of  this  implementation  with  the  implementation  of  Mod  I  under  UNIX  in  C  re¬ 
vealed  comparable  storage  efficiency  and  superior  access  efficiency.  Specifically,  the  bulletin 
and  arrival  print  functions,  running  on  the  same  day's  worth  of  data,  with  identical  environments 
(same  hardware,  single-user  mode),  ran  approximately  30-percent  faster  with  MUMPS  than 
with  UNIX.  In  terms  of  storage,  the  MUMPS  file  handler  required  approximately  65  kbytes  for 
the  sample  database,  compared  with  72  kbytes  on  the  UNIX  system.  The  MUMPS  figure  also 
includes  ancillary  inverted  files  created  for  the  query  capability  which  did  not  exist  in  the  UNIX 
file  structure. 

The  interpretive  character  of  the  MUMPS  language  makes  it  likely  that,  given  equal  levels 
of  skill  in  the  two  environments,  it  would  be  easier  and  faster  for  developing  the  parametric 
database.  The  main  issue  is  the  network  interface.  Because  the  available  implementations  of 
MUMPS  for  the  PDP-11  are  stand-alone  systems,  they  would  have  to  be  modified  to  handle  the 
network  interface  as  developed  by  Sytek.  Neither  of  the  two  vendors  (InterSystems  and  DEC)  is 
interested  in  performing  this  work.  For  this  reason,  it  is  apparent  that,  although  MUMPS  is 
an  attractive  environment 'for-the  parametric  database,  it  would  not  be  possible  to  use  it  in  the 
networked  environment  planned  for  the  SDC. 

H.  Lison 


G.  REMOTE  SEISMIC  TERMINAL 

The  Remote  Seismic  Terminal  (RST)  is  a  microcomputer-based  display  and  communication 
system  designed  to  support  the  analyst  faced  with  the  processing  of  data  from  a  seismic  station, 
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the  communication  of  seismic  data  to  the  SDC,  and  the  receipt  of  the  SDC  bulletin  at  a  remote 
location.  It  is  an  alternative  to  the  manual  processing  of  seismic  data  and  will  provide 
computer-assisted  communication  With  the  SDC.  In  its  basic  form,  the  RSY  will  provide  sup¬ 
port  to  the  analysis  of  digital  waveform  data  and  will  partially  automate  the  measurement  and 
extraction  of  seismic  parameters.  In  expanded  versions,  the  RST  can  digitize  waveform  data 
and  provide  a  variety  of  tools  for  the  analysis  of  digital  data.  An  overview  of  the  design  of  the 
RST  and  a  discussion  of  its  capabilities  are  given. 

The  RST  is  a  modular  microcomputer  and  peripheral  system  plus  a  modular  suite  of  soft¬ 
ware  designed  to  provide  communications  "dth  the  SDC,  waveform  display  capability  and  support 
for  creation  of  International  Exchange  of  Seismic  Data  (IESD)  Level  I  (parametric  data)  data 
reports. 

The  RST  hardware  consists  of  a  chassis  containing  a  microcomputer  and  associated  memory 
and  peripheral  boards,  an  alphanumeric  terminal,  a  graphics  terminal,  a  floppy  disk,  and  a 
modem  for  communications.  The  main  chassis  provides  power  supplies  and  the  means  for  in¬ 
terconnection  of  the  other  boards  via  the  standard  S-100  bus. 

In  the  main  chassis  are  four  boards  plus  the  memory  boards.  The  first  is  a  California 
Computer  Systems  (CCS)  CPU  board,  which  is  a  Z80  CPU.  Second,  the  Morrow  disk  controller 
board  operates  the  floppy  disks.  Third,  the  Scion  graphics  display  is  controlled  by  a  separate 
board  which  has  its  own  Z80  and  memory  connected  to  an  internal  bus.  Communication  with  the 
system  is  via  I/O  ports,  and  the  board  produces  video  for  driving  the  monitor.  The  final  main¬ 
frame  board  is  a  Morrow  Designs  "Switchboard"  which  incorporates  two  serial  and  four  parallel 
ports  plus  a  status  port. 

There  are  four  memory  boards  in  the  RST.  The  Z80  can  only  address  64K  of  memory  but 
the  S-100  bus  allows  the  switching  of  multiple,  independent  memory  banks  into  the  Z80  address 
space.  A  16K  static  RAM  is  used  for  the  stack  and  for  communication  between  programs  in  the 
other  memory  banks.  Next,  there  are  three  CCS  64K  dynamic  RAMs,  assigned  to  banks  1,  2, 
and  3,  respectively.  On  these  boards,  the  16  kbytes  of  each  are  disabled,  leaving  48K  of  each 
board  active.  The  bank  1  board  is  configured  to  become  active  on  power  on/reset,  and  thus  in 
"normal"  mode  the  entire  collection  looks  like  64K  of  contiguous  memory.  By  exercising  the 
bank  switching  logic,  the  48  K  bank  can  be  replaced  by  one  of  two  other  completely  independent 
sets. 

The  rest  of  the  hardware  consists  of  independently  packaged  units  interconnected  with  cables. 
The  dual-disk  drive  is  a  Morrow  Designs  Discus  2+2D  and  is  capable  of  reading  both  single-  and 
double-density  diskettes.  The  alphanumeric  terminal  is  a  standard  Ilazeltine  1421,  which  com¬ 
municates  with  the  system  through  a  serial  I/O  port.  A  Racal  Vadic  modem  is  used.  It  is  ca¬ 
pable  of  300-  or  1200-baud  operation,  the  latter  with  the  Bell  212A  (industry  standard)  protocol 
or  with  the  older  Vadic  protocol. 

The  software  suite  of  the  RST  includes  the  CP/M  (CP/M  is  a  registered  trademark  of  Dig¬ 
ital  Research)  operating  system  with  its  file  management  system,  editor,  compilers,  etc.  This 
supports  the  preparation  of  Level  I  reports  and  the  printing  of  bulletins.  It  also  supports  the 
generation  of  analysis  programs  and  the  analysis  of  data.  The  communications  function  is  sup¬ 
ported  by  specialized  programs  which  observe  a  communications  discipline  and  allow  the  trans¬ 
fer  of  files  of  daia  to  and  from  the  Data  Center.  Finally,  the  analysis  of  waveform  data  is  sup¬ 
ported  by  a  display  program  which  allows  the  display  of  and  measurement  from  waveform  data 
stored  on  the  local  floppy  disks. 
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The  RST  has  three  different  "operating  modes."  First,  it  is  or  can  be  a  normal  CP/M-baaed 
microcomputer  with  all  of  the  CP/M  features  such  as  editor,  assembler,  debugger,  peripheral 
interchange  processor,  and  so  on;  extra  features  are  a  BASIC  compiler  and  a  complete  relo¬ 
catable  macro  assembler  package  (Microsoft  M80), 

Next,  the  RST  can  be  configured  to  be  an  intelligent  terminal  for  use  in  communicating  with 
another  device  via  a  modem  supplied.  This  operating  mode  is  entered  from  the  CP/M  mode  by 
activating  an  appropriate  software  package. 

Finally,  the  RST  has  a  graphics  display  moSe,  again  entered  from  the  CP/M  mode  via  a 
call  to  DISP.  D1SP  is  a  self-contained  software  package  which  prompts  the  operator  to  perform 
various  waveform  display  and  analysis  operations,  DISP  takes  commands  from  the  terminal 
keyboard  and  allows  the  use  of  a  graphics  cursor  controlled  by  the  operator  to  extract  amplitude 
and  period  measurements  from  the  waveforms. 

A  number  of  different  enhancements  have  been  suggested  by  potential  users  who  have  had 
an  opportunity  to  see  the  RST  demonstrated  and  to  try  using  it  themselves.  They  group  into 
five  categories;  arithmetic  enhancement  through  the  addition  of  a  Pascal  processor;  addition 
of  a  hard  disk  for  increased  storage  capacity  and  access  speed;  software  enhancements  including 
hard  copy  of  graphics  display,  display  of  maps,  and  improved  communications  protocol;  addition 
of  a  plotting  board  for  digitization  of  analog  seismic  records;  and  real-time  operation  including 
direct  digital  data  acquisition  from  a  seismometer.  These  enhancements  are  being  planned  and 
should  be  accomplished  in  the  next  reporting  period. 

Considerable  interest  was  expressed  in  the  extent  to  which  the  demons i-ration  RST  config¬ 
uration  could  be  "cut  down."  In  Table  1-2,  five  systems  of  increasing  complexity  are  given  with 
their  costs.  Three  columns  of  cost  data  are  included;  the  first  column  is  the  incremental  cost 
of  adding  the  capability  described;  the  second  column  is  the  accumulation  of  these  increments, 
and  represents  the  cost  of  a  system  at  this  level  in  a  form  that  allows  future  expansion;  the 
third  column  shows  the  cost  of  an  all-in-one  system  at  the  desired  level  -  such  a  system  will 
be  cheaper  but  will  have  no  growth  capability. 


TABLE  1-2 

LEVELS  OF  COMP  LEX  fY  AND  COSTS 

Level  of  Capability 

Incremental 

Cost 

Total 

Cost 

Minimum  Cost 
Alternative 

Video  Terminal  +  Modem  +  Printer 

3,000 

3,000 

2,000 

+  Disk  +  CPU  +  Compiler 

4,500 

7,500 

5,000 

+  Graphics  Display  and  Cursor* 

2,500 

10,000 

(5,000)* 

+  Plotter  +  Hard  Disk§ 

5,000 

15,000 

(9,000)* 

+  Seismometer  Hookup 

21,000 

36,000 

N/A 

t  Model  demonstrated  at  the  Group  of  Scientific  Experts  Meeting. 
1  Limited  resolution  graphics,  limited  arithmetic  capability. 
§Currenl  level  of  capability. 
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The  first  system  allows  simple  communication  of  text,  it  has  no  computing  or  graphics 
capability.  An  even  smaller  system  is  possible  if  one  dispenses  with  the  lineprinter,  but  nobody 
was  willing  to  do  without  hard  copy.  The  next  increment  adds  a  computing  capability  (hardware 
and  software)  and  produces  a  viable  small  computer.  Addition  of  high -re solution  graphics  comes 
next;  the  (5,000)  in  the  third  column  is  intended  to  indicate  that  generally  the  preceding  system 
will  have  some  graphics  capability,  perhaps  100  X  280  or  better.  This  is  the  level  of  the  present 
RST.  The  next  two  increments  provide,  respectively,  a  realistic  database  capacity  with  analog 
data  reduction  capability,  and  a  true  seismic  station. 

R.  M.  Moroney 
A.  G.  Gann 


H.  VAX  UNIX  INSTALLATION 

In  December  1980  we  installed  the  Fourth  Berkeley  Release  of  Software  for  the  VAX*  on 
our  VAX  11/780  after  experimenting  for  a  short  time  with  VMS^  and  an  earlier  release  of  the 
Berkeley  material.  The  Fourth  Berkeley  Release  is  distributed  by  the  Computer  Systems  Re¬ 
search  Group  of  the  University  of  California  at  Berkeley,  and  consists  of  modifications  and  ex¬ 
tensions  to  Western  Electric's  UNIX/32V^  operating  system.  The  Fourth  Berkeley  Release  is 
largely  compatible  with  the  Version  7  UNirl  operating  system  which  is  currently  running  on 
our  PDP-il/50t  and  PDP-11/44. 

Despite  somewhat  lower  performance  in  particular  areas,  UNIX  was  selected  in  preference 
to  VMS,  the  operating  system  supplied  for  the  VAX  by  Digital  Equipment  Corporation  (DEC), 
because  it  offers  greater  ease  of  modification  for  special  functions,  a  more  complete  develop¬ 
ment  environment,  and  compatibility  with  our  existing  systems.  UC  Berkeley  has  been  funded 
by  DARPA  to  develop  a  standard  VAX  UNIX  system  for  Government  contractors,  and  to  enhance 
the  performance  of  UNIX  for  large  applications.  We  anticipate  that  the  system  to  be  produced 
from  this  project  will  make  UNIX  even  mere  desirable. 

The  only  substantial  obstacle  encountered  in  installing  this  system  was  a  lack  of  software 
support  for  our  UNIBUS^  interface  magnetic  tape  controller.  This  was  overcome  by  obtaining, 
from  UC  Berkeley,  an  experimental  software  driver  which  we  are  using  successfully. 

Initial  execution  speed  of  FORTRAN  code  compiled  under  UNIX  was  disappointing  when 
compared  with  FORTRAN  compiled  and  run  under  VMS,  Installation  of  a  hardware  floating¬ 
point  accelerator,  however,  resulted  in  a  dramatic  improvement  in  the  speed  of  UNIX  FORTRAN 
relative  to  that  of  DEC'S.  We  are  currently  awaiting  software  from  UC  Berkeley  '"hich  should 
improve  the  performance  of  UNIX  FORTRAN  even  further. 

D.  A.  Bach 


t  VAX,  VMS,  PDP,  and  UNIBUS  are  trademarks  of  Digital  Equipment  Corporation, 
t  UNIX  and  UNIX/32V  are  trademarks  of  Bell  Laboratories. 
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II.  SEISMIC  PROCESSING  TECHNIQUES 

A.  SEISMIC  ALGORITHM  DEVELOPMENT:  OVERVIEW 

Substantial  progress  has  been  made  in  the  development  of  the  seismic  algorithms  required 
for  Data  Center  operation.  Throughout  this  task,  the  emphasis  is  upon  replacing  the  functions 
currently  performed  by  the  human  analyst  by  automatic  algorithms,  while  retaining  the  capabil¬ 
ity  to  permit  analyst  intervention  for  the  revision  and  addition  of  the  automatically  determined 
parameters  at  several  points  during  the  data  flow  from  continuous  digital  data  to  final-event 
bulletin. 

We  have  previously  been  hampered,  particularly  with  respect  to  detection,  by  the  lack  of 
readily  available  data  for  testing  purposes.  This  problem  is  now  resolved  by  data  from  the 
International  Data  Collection  Experiment  (IDCE),  described  below,  which  provides  large  amounts 
of  data  suitable  for  testing  not  only  the  seismic  algorithms  but  also  many  other  aspects  of  SDC 
development.  Several  changes  have  been  made  to  the  detection  algorithm  which  is  described  in 
some  detail,  and  a  description  of  its  performance  upon  one  day  of  IDCE  data  is  given.  A  pro¬ 
cess  to  revise  the  signal  onset  times  provided  by  the  detector  and,  additionally,  to  identify 
source  character  and  measure  magnitude-related  parameters  is  described  and  evaluated.  A 
multiple-filter  technique  to  time  and  ii'entify  body  waves  is  also  outlined  and  applied  to  local 
data.  We  emphasize  that  any  application  of  the  above  algorithms  for  short-period  data  may  well 
require  substantial  modification  when  the  broader-band  high-frequency  RSTN  data  become  avail¬ 
able.  A  proposed  method  for  "sym  rtic"  analysis  of  seismograms,  designed  to  recognize  and 
use  the  features  of  seismic  data  employed  by  the  analyst,-  is  described  and  its  possible  applica¬ 
tions  are  discussed. 

The  low  noise  characteristics  of  the  SRO/ASRO  long-period  data  imply  that  surface-wave 
detection  capability  may  be  superior  to  that  given  by  short-period  body  waves,  particularly  for 
regions  with  poor  station  coverage  such  as  the  southern  hemisphere.  Techniques  for  the  detec¬ 
tion,  azimuth  determination,  and  association  of  Rayleigh  waves  to  provide  a  location  based  upon 
long-period  data  only  are  described  and  tested  using  IDCE  data. 

R.  G.  North 


B.  IDCE  DATABASE 

From  i-15  October  i960,  several  countries  cooperated  in  the  International  Data  Collection 
Experiment  proposed  by  Sweden  to  provide  experience  in 

(1 )  the  mechanics  of  an  international  exchange  of  seismic  data, 

(2)  the  extraction  of  Level  I  parameter  data,  and 

(3)  the  handling  of  Level  I  and  il  data  at  Data  Centers. 

Parameter  and,  where  available,  waveform  data  were  recorded  at  many  stations  worldwide 
and  sent  to  Sweden  where  the  data  will  be  collated  and  a  bulletin  prepared.  The  data  from  this 
experiment  will  be  invaluable  for  testing  all  functions  of  the  SDC,  and  will  be  the  only  data  avail¬ 
able  for  this  purpose  until  the  RSTN  starts  operation. 

The  U.S.  contribution  to  the  IDCE  consisted  of  digital  data  recorded  at  2  5  stations  whose  lo¬ 
cations  are  shown  in  Fig.  Il-i.  This  network  includes  the  Global  Digital  Seismograph  Network 
(GDSN )  comprising  the  SRO,  ASRO  and  DW;WSSN  sites  operated  by  theUSCS.  Additional  stations 
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include;  USGS  Seismic  Network  Bites  A  OK,  BDW,  PM  11,  and  SMYj  SDCS  sites  P'-WY  and 
HN-IYK;  and  the  prototype  HSTN  station  at  CPO.  Twenty  of  the  stations  (9  of  which  are 
3-coinponent)  recorded  continuous  short-period  datu,  while  the  remaining  5  used  an  on-site  de¬ 
tector  to  select  event  waveforms,  All  stations  recorded  continuous  3-component  long-period 
data.  While  the  inevitable  data  transmission  dropouts  and  occasional  station  failures  occurred, 
to  our  knowledge  this  is  the  largest  digital  database  yet  collected  from  a  global  network.  Fur¬ 
ther  data  contributed  by  other  nations  will  eventually  become  available  and  augment  thiB  database. 

In  addition  to  the  Level  II  waveform  data,  Level  1  parameters  were  also  supplied  to  Sweden 
by  the  U.S.  These  consist  of  analyst  m<  surements  for  most  of  the  2  5  digital  stations,  made 
from  the  analog  seismograms  which  are  recorded  at  each  site.  For  each  time  pick,  the  maxi¬ 
mum  amplitude  and  dominant  period  were  measured  by  L’SGS  analysts,  and  a  special  effort  was 
made  to  Identify  secondary  phases  and  corresponding  arrivals  on  the  long-period  components 
and,  where  available,  the  horizontal  component  short-period  data.  These  analyst  picks  will  be 
invaluable  in  evaluating  automatic  detection  and  signal  measurement  algorithms  for  the  SDC. 

On  2  and  4  October  the  time,  amplitude,  and  period  measurements  were  repeated  with  the  digi¬ 
tal  data  in  order  to  compare  these  with  analog  determinations. 

The  above  is  a  brief  summary  of  a  detailed  report1  prepared  by  the  USGS  for  presentation 
at  the  February  1 981  meeting  in  Geneva  of  the  Group  of  Scientific  Experts  of  the  Committee  on 
Disarmament.  We  have  received  from  SDAC  copies  of  the  digital  data  tapes  and,  from  USGS/ 
Golden,  the  parameter  database.  A  major  advantage  of  the  digital  data  is  their  standardization 
to  the  SRO  tape  format.  While  some  minor  modifications  to  this  format  were  unavoidably  re¬ 
quired  to  accommodate  the  range  of  data  types,  we  have,  with  gratefully  acknowledged  assis¬ 
tance  from  both  USGS/Golden  and  SDAC,  been  hie  to  read  digital  data  from  all  the  stations 
with  small  changes  to  our  program  to  read  the  SRO  network  day  tapes.  This  program  converts 
the  data  into  the  waveform  database  structure  which  we  have  developed,  and  for  which  we  have 
written  much  special-purpose  software  for  the  display  and  analysis  of  digital  data.  We  are  thus 
able  to  immediately  start  data  analysis  and  testing  of  the  SDC  seismic  algorithms  using  this 
large  database. 

The  digital  data  consume  over  100  tapes  at  1600  bpi,  and  data  volumes  are  ~70  million 
samples/day,  or  -280  Mbytes  when  stored  in  the  floating-point  representation  used  by  the  wave¬ 
form  database  structure.  This  exceeds  Ly  nearly  2  5  percent  the  data  rates  anticipated  from  the 
5-station  RSTN.  and  its  processing  will  enable  reliable  estimates  to  be  made  of  the  computa¬ 
tional  requirements  and  special  hardware  needs,  such  as  array  processors,  •■•ecessary  for  an 
operational  version  of  the  SDC. 

We  have  prepared  sample  databases  in  the  waveform  database  format  for  the  2  and  4  Octo¬ 
ber  data.  At  present,  we  do  not  have  the  disk  storage  capacity  to  maintain  one  day  of  continu¬ 
ous  data  from  all  stations  in  on-'ine  storage,  though  wo  anticipate  possessing  the  capability  to 
store  several-day's  worth  on-line  in  the  very  near  future.  The  contiruous  data  have  thus  been 
initially  stored  on  UNIX  Tar  (Tape  Archive)  tapes,  from  which  retrieval  is  much  more  rapid 
than  from  the  original  dat.  tapes.  We  have  run  our  detection  algorithm  on  all  the  continuous 
data  for  the  two  days  selected  and  intend  to  maintain  the  detected  waveform  segments,  together 
with  the  Level  I  parameters,  on-line  indefinitely.  This  will  constitute  the  major  data  source 
for  testing  analyst  procedures  such  as  detection  review  and  revision  of  automatically  measured 
parameter  determinations  upon  the  HAS. 
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The  data  quulity  is  generally  good,  particularly  at  the  GDSN  sites.  Stations  whose  data 
were  telemetered  ."om  the  site,  such  as  the  Alaskan  IJSGS  stations  and  the  SDCS  sites,  Include 
numerous  data  dropouts  and  glitches.  We  have  encountered  some  timing  problems  at  CPO  and 
some  of  the  telemetered  stations.  At  BDW  (Boulder,  Wyoming)  there  is  apparent  data  quanti¬ 
zation  due  presumably  to  partial  failure  of  the  A/D  converter.  While  these  assorted  features 
are  undesirable,  they  are  representative  of  the  problems  with  which  the  SDC,  and  the  seismic 
algorithms  in  particular,  must  be  able  to  cope.  Later  in  this  report  we  describe  some  pre¬ 
liminary  results  of  analysis  of  these  data,  particularly  in  the  ureas  of  detection  and  parameter 


measurement. 


R.G.  North 


C.  DETECTION  ALGORITHM 

The  detection  algorithm  which  we  have  implemented  and  described  very  briefly  in  previous 
SATSs^*^  is  basically  the  Z-statistic  detector^  Lith  a  few  ad., -'tons.  We  have  recently  gener¬ 
alized  die  detector  to  permit  rapid  change  of  parameters,  added  »  glitch  remover,  and  Included 
the  recursive  bandpass  prefilter  in  the  detection  program  itself.  While  the  Z-statistic  detector 
has  been  described  and  evaluated  in  some  detail  by  its  original  developers,^  we  present  here  a 
summary  of  its  implementation  by  this  group, 

A  sample  input  parameter  file  is  given  below: 


0.50 

low  corner  freq 

(hz) 

(1) 

4.00 

high  corner  freq 

(hz) 

(2) 

0.00 

mu 

<  log  <  x*  #  2>i 

(3) 

0.00 

sigma 

(var)*  *  .5 

(4) 

60,00 

LTA  window 

(sec) 

(5) 

3.00 

STA  window 

(sec) 

(6) 

1.50 

STA  increment 

(sec) 

'7 ) 

20.00 

STA  to  LTA  lag 

(sec) 

(8) 

50.00 

coda  reset 

(sec) 

(9) 

3.75 

detection  threshold 

(10) 

1 

minimum  detection  votes 

(11) 

yes 

freeze  LTA  during  detections? 

(12) 

yes 

de-glitch  ? 

(13) 

szr 

component 

(14) 

In  the  above,  (1 )  and  (2)  are  the  3-dB  low-  and  high-frequency  cutoffs  of  the  prefilter  band¬ 
pass.  The  present  filter  is  a  3-pole  Butterworth;  faster  roll-off  than  this  seems  desirable  and 
will  be  implemented  shortly.  Lines  (3)  and  (4)  are  initial  values  of  p  and  <j  for  log  (STA),  where 
STA  is  the  short-term  power  over  a  time  window  of  given  length  (here  3.0  s).  The  statistics  p 
and  a  are  here  set  to  zero  for  a  "cold  start,"  but  in  routine  operation  use  of  previous  values  for 
the  same  site  would  eliminate  the  warmup  time  which  this  choice  requires.  Line  (5)  gives  the 
long-term  averaging  time  constant  which  determines  how  rapidly  new  values  of  STA  are  incor¬ 
porated  in  the  update  of  p  and  a.  — - - — . . 

The  STA  need  not  be  updated  nor  its  Z  tested  against  the  threshold  for  every  sample  point, 
here  [line  (7)]  it  is  shifted  by  1.50  s  or  half  its  window  length  each  time.  It  could  be  shifted  by 
the  entire  window  length,  but  this  raises  the  possibility  that  a  small  signal  may  be  split  between 
successive  windows,  neither  of  which  has  sufficient  power  to  trigger  the  detector.  The  STA  to 
LTA  lag  [line  (8)|  was  introduced  to  assist  detection  of  emergent  arrivals,  which  would  slowly 
pollute  p  a  id  e  so  that  the  threshold  would  never  be  exceeded.  The  lag  here  (20  s)  is  the  time 
between  the  STA  being  tested  and  the  STA  to  be  used  in  LTA  update.  Lags  of  10  to  20  s  have 
been  found  to  optimize  performance  for  most  emergent  arrivals.  The  coda  reset  time  [line  (9)J 
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TABLE  11-1 

SAMPLE  OUTPUT  OF  DETECTION  LOG 


0 

I 

DETECTOR  LOG  onto4 

1 

1 

I 

low  comer  freq  =  ,50  hz 

I 

1 

high  corner  freq  =  4,00  hz 

0 

1 

mu  =  .00 

0 

1 

sigme  =  .00 

0 

# 

LTA  window 

=  60.00  sec 

0 

1 

STA  window 

«  3 .00  sec  ■ 

0 

I 

STA  Increment  “  1 .50  see 

0 

I 

STA  to  LTA  leg  -  20.00  toe 

0 

# 

coda  reset  = 

50.00  sec 

0 

1 

detection  threshold  >=  3.75 

0 

1 

minimum  det#ct*on  votes  ■ 

1 

0 

If 

freeze  ■  yes 

0 

0 

de-glitch  = 

yes 

0 

1 

0 

component  = 

izr 

1 

1 

# 

Gram  1 

23:59:59.960 

0 

1  10/  4/80 

1:21:46.410 

start 

DETECTOR 

1  10/  4/80 

1:22:38.910 

end 

Zmax  *  4.0 

1  IQ/  1/80 

2:28:14.910 

start 

DETECTOR 

1  10/  ',/80 

2:29:  8.910 

end 

Zmax  *  4.7 

1  IQ/  4/80 

3:42:20.910 

stort 

DETECTOR 

1  10/  4/80 

3:43:20.910 

end 

Zmax  =  4.3 

1  10/  4/80 

4:57:20.910 

start  DETECTOR 

1  10/  4/80 

5:  0:23.910 

end 

Zmax  -  9.4 

1  10/  4/80 

12:22:11.911 

start  DETECTOR 

1  10/  4/80 

12:23:  4.411 

end 

Zmax  «  4.1 

1  10/  4/80 

15:13:11.911 

start 

DETECTOR 

1  IQ/  4/80 

15:20:44.911 

end 

Zmax  =  32.4 

1  IQ/  4/80 

16:  6:47.91 1 

stort 

DETECTOR 

1  10/  4/80 

16:10:50.911 

end 

Zmax  =  17.4 

1  IQ/  4/80 

16:29:17.911 

start 

DETECTOR 

1  IQ/  4/80 

16:30-19.41 1 

end 

Zmax  =  7.9 

1  10/4/80  16:38:29.91  T 

start  DETECTOR 

1  10/  4/80 

16:39:23.911 

end 

Zmax  =  4,6 

1  10/  4/80 

19:  2:41.911 

start 

DETECTOR 

1  10/  4/80 

19:  3:34.411 

end 

Zmax  =  4.5 

1  10/  4/80 

19:  6:25.411 

start 

DETECTOR 

1  10/  4/80 

19:  7:37.411 

end 

Zmax  =  12.4 

1  10/  4/80  20:1 1:22 .411 

start 

DETECTOR 

1  10/  4/80  20:12:14.911 

end 

Zmax  =  3 .8 

1  10/  4/80  23:54:50.911 

start 

DETECTOR 

1  10/  4/80  23:55:43.411 

end 

Zmax  =  3.8 

# 

EOF 

0:  0:35.861 

0 

# 

mu  =  18.15 

0 

# 

sigma  = 

35 

0 

is  the  minimum  amount  of  time  allowed  between  detections,  and  inhibits  multiple  detections 
within  a  coda.  For  windowing  purposes,  this  time  is  considered  part  of  the  detection  interval. 

The  minimum  detection  votes  [line  (14  )J  may  be  used  in  an  attempt  to  prevent  detection  of 
glitches  and  anomalously  short-lived  arrivals.  If  set  to  3,  for  example,  then  it  requires  that 
the  threshold  be  exceeded  at  least  3  times  in  succession  before  a  detection  is  declared.  We 
have  found  its  usefulness  to  be  low.  The  t.TA  (i.e.,  p.a)  may  be  "frozen*  (line  (12))  when  a  de¬ 
tection  is  declared:  this  prevents  their  contamination  by  signal  characteristics  and  reduces  the 
number  of  multiple  detections  during  a  given  arrival.  On  the  negative  side,  freezing  may  in¬ 
hibit  detection  of  signals  from  another  event  during  u  detection  interval,  but  this  function  may 
be  best  performed  by  the  analyst  during  detection  review. 

The  de-glitch  option  [line  (13))  permits  application  of  a  simple  technique  to  remove  large 
single-point  excursions.  Prior  to  prefiltering,  the  mean  is  removed  and  the  absolute  value  of 
each  point  is  tested  against  that  of  its  immediate  neighbors.  If  the  former  exceeds  the  latter 
by  more  than  a  factor  of  10,  it  is  deemed  a  glitch  and  set  to  the  average  value  of  its  neighbors. 
We  plan  to  test  some  ideas  for  a  more  sophisticated  post-detection  multipoint  de-glitch  routine, 
and  candidates  include  checks  for  excessively  one-sided  arrivals  and  large  ratios  of  the  max¬ 
imum  Z  during  detection  to  the  detector  on  time. 

The  ability  to  set  all  these  relevant  parameters  permits  rapid  implementation  of  changes 
to  improve  performance.  A  sample  output  detection  log  for  the  detector  running  with  the  param¬ 
eters  given  above  for  24  h  of  short-period  vertical  data  at  station  ANTO  is  given  in  Table  II— 1 . 
The  chosen  detector  parameters  are  given,  together  with  the  start  and  end  times  of  the  detec¬ 
tions.  For  each  detection,  the  maximum  value  of  Z  within  the  detection  interval  is  given.  This 
permits  rapid  identification  of  the  signals  which  would  no  longer  be  detected  if  the  threshold 
were  raised  by  a  specified  amount.  This  sample  run,  including  glitch  removal  and  prefiltering, 
took  only  ~10  min.  on  the  VAX  computer  and,  thus,  the  effects  of  parameter  changes  can  be 
rapidly  investigated. 

A  special-purpose  windowing  program  has  been  written  to  ren^  the  deit-i^on  log  and  extract 
windows  corresponding  to  the  detection  intervals.  A  selectable  additional  amount  'f  time  is 
added  before  and  after  the  start  and  end  times  indicated  when  the  window  is  extracted.  This  is 
currently  set  at  60  s.  Waveforms  corresponding  to  the  above  intervals  are  shown  in  Fig.  II-2. 

An  evaluation  of  the  detection  algorithm  using  IDCE  data  and  a  discussion  of  the  scope  for 
improvements  which  this  experiment  has  revealed  are  presented  in  Sec.  D  below. 

M.  A.  Tiberio 


D.  DETECTOR  EVALUATION  WITH  IDCE  DATA 

We  have  applied  the  detection  algorithm  described  above  to  continuous  short -period  vertical 
component  data  at  1 9  stations  provided  by  the  IDCE  for  4  October  1980.  A  valuable  check  upon 
its  performance  is  provided  by  the  analyst  picks  made  for  this  period  and  the  NEIS  bulletin  which 
gives  16  events  located  for  this  day. 

Table  II— 2  lists  the  detection  parameters  used  -  these  are  described  in  more  detail  in  Sec.  C 
above.  We  decided  to  apply  the  same  detector  parameters  to  all  stations  for  consistency,  though 
in  practice  different  prefilters  and  short-term  power  windows,  selected  on  the  basis  of  experi¬ 
ence,  will  provide  optimum  performance  at  each  sitj.  We  initially  selected  the  band  1  to  4  Hz 
for  the  prefilter  bandpass,  but  discovered  that  this  choice  resulted  in  a  failure  to  detect  some 
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TABLE  11-2 

DETECTOR  PARAMETERS  USED  FOR  ALL  STATIONS 

0.50 

low  corner  freq  (hi) 

4.00 

high  corner  freq  (hi) 

0.00 

mu  <lnO**2» 

0.00 

slgmo  (vor)**.5 

60.00 

LTA  time  constant  (see) 

3.00 

STA  window  (sec) 

1.30 

STA  Increment  (see) 

20.00 

STA  to  LTA  lag  (sec) 

50.00 

code  reset  (sec) 

t  3.85 

detection  threshold 

1 

minimum  detection  votes 

y« 

freeze  LTA  during  detections? 

y«* 

deglltch? 

sir 

component 

*  Adjusted  for  optimum  performance  at  each  station. 

signals  (particularly  PKP)  from  the  largest  event  on  this  day.  This  failure  may  not  be  particu¬ 
larly  ser,ous  since  the  number  of  P  signals  detected  from  this  event  would  still  be  adequate  for 
location  purposes,  but  we  decided  to  change  the  prefilter  bandpass  to  0.5  to  4  Hz.  This,  in  its 
turn,  resulted  in  a  failure  to  detect  some  smaller  high-frequency  signals  from  local  events  with¬ 
out  setting  the  threshold  for  detection  so  low  that  unacceptably  large  false-alarm  rates  resulted. 
We  believe  that  this  may  indicate  a  need  to  run  the  detector  several  times  upon  the  same  data, 
with  perhaps  as  many  as  three  different  bandpass  prefilters;  and  while  this  can  easily  be  carried 
out  with  the  current  detection  algorithm,  the  logic  for  combining  detections  in  a  number  of  such 
bands  into  a  final  detection  log  is  not  straightforward  and  will  be  investigated  in  the  near  future. 

Detector  performance  proved  to  be  relatively  insensitive  to  short-term  power  window  lengths 
in  the  range  1  to  3  s  and  also  to  long-term  averaging  times  in  the  range  20  to  60  s.  The  detec¬ 
tion  threshold  was  adjusted  for  each  station  to  provide  optimum  detection  of  the  analyst-picked 
signals  combined  with  a  moderate  false-alarm  rate  which  was  chosen  to  result  in  -30  detections/ 
day  at  most  stations.  Fewer  detections  were  made  at  the  high-noise  stations  CTAO.  GllMO,  and 
KQNO.  A  very  large  number  of  detections  were  made  at  CPO,  but  these  are,  without  exception, 
signals  from  presumed  mine  and  quarry  blasts  at  distances  of  up  to  300  km.  No  teleseisms 
were  detected  at  this  station  due  to  a  combination  of  this  almost-continuous  local  activity  and 
extremely  high  bad  ground  noise  levels. 

The  results  of  the  detection  processing  and  its  comparison  with  analyst  performance  are 
summarized  in  Table  II-3.  For  each  station  the  number  of  detections  and  the  percentage  of  time 
the  detector  was  on  during  the  24-h  interval  are  given.  This  exceeds  4  percent  only  at  CPO. 

The  total  percentage  of  data  saved  is  larger,  since  we  decided  to  window  the  detections  for  an¬ 
alyst  review  with  an  additional  60  s  before  and  after  the  detector  start  and  end  times,  respec¬ 
tively.  This  is  necessary  to  ensure  that  P  is  included  for  local  arrivals  since  the  detector 
often  triggered  on  the  stronger  S  phase  at  these  distances,  and  the  extension  of  the  window 
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TABLE  11-3 


DETECTION  STATISTICS  AND  COMPARISON  WITH  ANALYST  PICKS 


Station 

Pick 

Detections 

Miss 

Total 

Detections 

Extra 

Real 

Folse 

Alarm 

Glitch 

Percent 

On 

Percent 

Save 

ANMO 

13 

12 

n 

25 

13 

8 

5 

0 

2.3 

5.5 

ANTO 

6 

6 

13 

7 

4 

3 

0 

1.7 

3.0 

BCAO 

9 

8 

n 

26 

18 

15 

3 

0 

3.1 

5.9 

BOCO 

2 

2 

13 

11 

1 

0 

1.2 

2.5 

GRFO 

6 

6 

n 

22 

16 

8 

0 

2.0 

4.3 

GUMO 

4 

3 

i 

24 

21 

5 

0 

2.2 

JASt 

1 

1 

1 

8 

7 

1 

0 

0.3 

NWAO 

1 

0 

i 

3 

3 

3 

0 

0.2 

ZOBO 

17 

14 

3 

31 

17 

10 

7 

0 

3.5 

59 

52 

B 

165 

113 

77 

36 

0 

1.8 

3.8 

ADK 

25 

3 

22 

31 

28 

1 

B 

20 

2.3 

5.5 

BDW 

18 

15 

3 

35 

20 

3 

11 

3.8 

7.5 

PMR 

21 

18 

3 

43 

26 

7 

I 

16 

3.7 

8.2 

SMY 

1 

1 

0 

33 

32 

10 

B 

15 

2.6 

6.1 

65 

37 

28 

142 

106 

21 

23 

62 

3.1 

6.8 

CPO 

2 

1 

B 

92 

91 

80 

11 

Bi 

8.6 

18.? 

CTAO 

- 

- 

- 

13 

13 

4 

Bl 

1.2 

2.1 

KONO 

- 

- 

2 

2 

0 

SI 

0.1 

0.3 

MAJO 

- 

- 

17 

17 

14 

3 

0 

1.8 

3.6 

32 

32 

25 

B 

0 

■ 

1.0 

2.0 

18 

18 

0 

46 

28 

3 

25 

0 

- 

5.6 

KAAO 

15 

10 

5 

30 

20 

15 

3 

2 

- 

6.1 

LON 

10 

9 

1 

22 

13 

13 

0 

0 

- 

5.1 

SNZO 

1 

0 

1 

3 

3 

0 

3 

0 

- 

2.6 

TATO 

3 

1 

2 

3 

2 

0 

2 

0 

- 

0.5 

47 

38 

9 

104 

66 

31 

33 

2 

D 

ffor  JAS,  continuous  data  were  recorded  only  for  the  first  10.2  h  and  the  statistics 
refer  to  this  portion  of  the  day  only. 


past  the  switch-off  time  of  the  detector  is  designed  to  include  the  pP  depth  phase  for  events  up 
to  400  km  in  depth.  While  it  might  be  desirable  to  increase  the  time  saved  past  the  end  of  the 
detection  to  include  all  conceivable  depth  phases  (240  s  required)  for  teleseisms,  and  Lg  for 
regional  events  (300  s  for  up  to  15°  distance),  this  would  greatly  increase  the  amount  of  data 
saved.  The  continuous  data  are  always  available  for  rewindowing  to  enable  the  analyst  to  search 
for  these  additional  later  arrivals. 

Table  II-3  also  gives  the  number  of  initial  (P)  arrivals  picked  by  the  analyst  for  each  station, 
and  the  number  of  these  detected  and  missed  by  the  detector.  The  number  of  extra  detections 
not  picked  by  the  analyst  are  also  given:  these  include  glitches  and  false  alarms,  but  many  are 
genuine  signals.  Detector  performance  was  excellent  at  most  of  the  SRO/ASRO  stations,  but 
less  acceptable  at  USGS  stations  ADK,  BDW,  and  PMR.  In  the  latter  case,  much  of  this  is  due 
to  the  excellent  analyst  performance  at  these  sites.  Many  of  the  signals  missed  by  the  detector 
can  be  captured  by  tailoring  the  prefilter  bandpass  to  each  individual  station.  Figure  II-3  shows 
examples  of  the  analyst-picked  signals  which  evaded  the  detector  with  the  standard  (0.5  to  4  Hz) 
prefilter  bandpass.  Analyst  picks  for  signals  5,  6,  and  10  are  PKP  arrivals  from  an  m^  5.5 
event  in  the  Fiji  region  at  noisy  sites.  At  NWAO  (6)  and  GUMO  (10),  the  arrivals  are  only  barely 
distinguishable  from  the  noise  and  the  onset  time  is  difficult  to  identify.  At  CPO  (5),  we  fail  to 
see  how  a  pick  could  have  been  made.  Local  arrivals  2,  3,  7,  and  9  can  easily  be  detected  by 
changing  the  prefilter  to  the  1-  to  4-Hz  bandpass,  with  scarcely  any  increase  in  false-alarm 
rate.  At  BDW  (4),  it  can  be  seen  that  the  data  are  quantized,  due  presumably  to  either  partial 
failure  or  incorrect  setting  of  the  A/D  converter,  and  we  are  surprised  that  under  these  cir¬ 
cumstances  we  have  been  able  to  detect  1  5  of  the  18  signals  picked  at  this  station.  The  appar¬ 
ently  disastrous  failure  to  detect  signal  8  at  ZOBO  can  be  explained  by  its  occurrence  only  100  s 
after  the  start  of  the  day,  before  the  noise  statistics  have  been  initialized.  We  shall  attempt  to 
correct  this  type  of  failure  due  to  "warm-up"  time.  The  worst  performance  is  that  for  ADK,  at 
which  only  3  of  the  2  5  picks  were  detected.  Of  the  31  detections  made,  20  were  small  glitches 
which  evaded  the  de-glitching  algorithm.  These  glitches  were  not  only  small  single-point  data 
excursions,  but  also  small  one-sided  pulses  which  were  clearly  not  of  seismic  origin.  Nearly 
all  the  signals  missed  were  very  high-frequency  local  events,  and  by  changing  the  prefilter 
bandpass  to  be  2  to  6  Hz,  20  of  the  25  picks  were  detected,  but  with  about  a  200-percent  increase 
in  false  alarms  which  were  again  mainly  glitches. 

On  the  brighter  side,  many  detections  were  made  which  did  not  correspond  to  analyst  time 
picks.  At  the  SRO/ASRO/DWWSSN  stations,  a  large  number  of  these  (—70  percent)  are,  in  our 
opinion,  genuine  seismic  arrivals  for  which  time  picks  can  be  made.  At  the  stations  ADK,  BDW, 
PMR,  and  SMY  nearly  all  were  either  false  alarms  or  glitches.  Figure  11-4  shows  some  ex¬ 
amples  of  detections  of  signals  which  were  not  picked  by  the  analyst.  For  the  stations  in  com¬ 
mon  between  this  figure  and  Fig.  II  —  3  (undetected  picks),  we  feel  that  the  unpicked  detections 
are  in  general  stronger  than  the  undetected  picks.  Four  of  the  signals  shown  in  Fig.  11-4  in  fact 
correspond  to  arrivals  predicted  by  the  PDE  event  list. 

The  last  five  stations  given  in  Table  II- 3  were  those  for  which  continuous  data  were  not  re¬ 
corded,  and  the  statistics  shown  refer  to  the  on-site  detector  at  these  sites.  The  performance 
of  this  detector  compared  with  the  analyst  was  exemplary  at  CIITO,  but  a  large  number  of  the 
unpicked  detections  were  false  alarms. 

To  give  some  indication  of  the  ability  to  locate  events  using  arrivals  from  an  automatic 
detector  only,  with  the  analyst  intervening  only  to  revise  the  onset  time  given  by  the  detector, 
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Table  II-4  shows  the  number  of  detections  made  for  each  of  the  16  events  given  in  the  PDE  for 
4  October.  Of  these  events,  5  were  recorded  only  at  local  stations  and  were  assigned  no  mag¬ 
nitude.  Of  the  remaining  11,  it  is  our  opinion  that  the  azimuthal  coverage  of  the  detections  made 
would  have  been  sufficient  for  a  prePminary  location  to  be  made  in  7  cases.  We  are,  to  say  the 
least,  puzzled  by  the  lack  of  detections  (and  analyst  picks)  for  the  m^  5.2  event  in  the  Philippines. 
The  only  arrival  corresponding  to  the  predicted  arrival  times  from  this  event  is  that  at  KAAO. 


TABLE  11-4 

EVENTS  GIVEN  IN  THE  NEIS  MONTHLY  EVENT  SUMMARY 
FOR  4  OCTOBER  1980 


Origin  Time 

Location 

"b 

Depth 

00:58:26.2 

Greece 

4.0 

10 

02:08:43.8 

Tonga 

4.7 

107 

03:11:35.8 

Columbia 

0 

167 

03:24:49.0 

Philippines 

5.2 

33 

04:10:31.0 

Chile 

4.7 

33 

04:37:34. 1 

Fiji 

5.5 

33 

06:09:03.3 

Italy 

0 

10 

06:32:13.7 

Sumatra 

4.6 

33 

12:56:59.6 

New  Hebrides 

4.8 

208 

15:12:03.8 

Greece 

4.9 

10 

16:05:40.0 

Turkey 

3.8 

10 

16:37:12.7 

India 

4.5 

33 

16:38:22.6 

California 

0 

4 

16:42:18.4 

California 

0 

5 

22:46:16.3 

California 

0 

5 

23:35:43.6 

Japan 

4.9 

366 

Stations  Detecting 


10  PMR? 


BDW,  PMR 

BOCO 

KAAO 

AMMO,  BCAO,  BDW,  ZOBO 

ADK,  ANMO,  ANTO,  BCAO,  CHTO, 
CTAO,  GRFO,  JAS,  LON,  MAJO, 
BDW,  PMR,  SMY,  ZOBO 


BCAO,  CHTO,  CTAO,  MAJO 

ANTO,  BCAO,  CHTO,  GRFO,  KAAO, 
KONO,  MAJO,  PMR,  SMY 

ANTO,  BCAO,  GRFO 

BCAO,  CHTO,  KAAO 

ANMO,  BDW 


ANMO,  CHTO,  GRFO,  GUMO,  KAAO, 
LON,  MAJO,  TATO,  ZOBO 


From  our  experiences  with  these  date,  we  have  gained  valuable  insights  into  methods  for 
improving  detector  performance.  Specific  changes  which  would  improve  detection  include: 


(1)  Detection  in  several  frequency  bands, 

(2)  Prefilter  optimized  to  station  noise  characteristics, 

(3)  Better  de-glitcher, 

(4)  Special  techniques  during  detector  "warm-up"  to  avoid  missing  signals 
during  this  interval. 

Additional  changes  which  we  plan  to  implement  shortly  include  3-component  detection  and 
methods  tailored  to  arrivals  of  longer  duration  such  as  local  S  and  regional  Lg. 


R.  G.  North 
M.A.  Tiberio 


E  E.  AUTOMATIC  TIMING  AND  ANALYSIS  OF  DIGITAL  DATA 
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For  each  signal  detected,  the  uetector  described  in  the  previous  report  outputs  a  windowed 
seismogram  containing  the  detected  signal  and  60  s  of  the  time  series  preceding  it,  and  a  marker 
indicating  the  detector  onset  time.  This  is  used  as  input  to  a  second  program  which  attempts  to 
refine  the  signal  onset  time  and  measure  other  seismic  parameters  from  the  signal. 

There  are  several  advantages  to  this  two-step  approach.  Since  the  detection  and  analysis 
steps  can  be  developed  independently,  the  software  is  less  complicated  and  easier  to  develop. 
Also,  each  step  can  have  different  design  goals.  The  detector  can  be  designed  to  be  simple  and 
efficient.  The  analyzer  can  be  designed  to  be  more  thorough.  Since  a  large  number  of  detected 
waveform  segments  can  be  kept  on-line,  the  development  of  the  analyzer  is  simplified. 

The  analysis  program  under  development  has  three  main  steps.  First,  it  uses  the  frequency 
content  of  the  signal  to  choose  an  appropriate  prefilter  bandpass.  Second,  an  arrival  time  pick¬ 
ing  algorithm  similar  to  that  of  Allen5  is  used  to  estimate  the  onset  time.  Finally,  the  onset 
time  is  refined,  the  amplitude  and  period  are  measured,  and  the  arrival  time  accuracy  is  esti¬ 
mated.  Each  of  these  steps  will  be  described  more  fully  in  turn. 

The  ratio  of  high-frequency  (3  to  8  Hz)  to  low-frequency  (0.3  to  1.0  Hz)  power  in  the  5  s  fol¬ 
lowing  the  detector's  trigger  time  is  used  to  filter  the  seismogram  and  decide  if  it  is  local  or 
telesuismic.  If  the  ratio  is  greater  than  0.1,  the  trace  is  high-pass  filtered  using  a  3-pole 
Butterworth  filter  with  the  low-frequency  3-dB  point  at  2  Hz.  Otherwise,  it  is  Butterworth 
bandpass  filtered  with  3-dB  points  at  0.667  and  3.0  Hz.  The  trace  is  only  forward  filtered. 
Forward  filtering  with  this  filter  causes  a  phase  delay  that  we  avoid  by  forward  and  backward 
filtering.  However,  the  sharper  phase-free  filtering  that  results  pushes  energy  from  impulsive 
arrivals  backward  in  time,  and  causes  the  program  to  pic  k  a  first  arrival  that  is  too  early. 

For  a  population  of  events  with  known  epicenters  recorded  at  MAJO,  North  and  Shields5 
found  events  with  a  ratio  of  0.1  or  greater  to  be  local  (<  500  km).  Events  with  a  ratio  of  0.03  or 
less  were  found  to  be  teleseisms.  Adjustments  to  this  test,  if  necessary  for  other  stations, 
will  be  incorporated  into  the  program  as  information  becomes  available. 

Once  the  data  have  been  filtered  and  the  event  roughly  categorized  by  distance,  a  picker 

5 

similar  to  that  of  Allen  is  used  to  estimate  the  signal  onset  time.  Allen  defines  the  E  function 
as 


E(t)  =  f(t)2  +  [C2  *  f'(t)]2 
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( 


where  f(t)  is  the  seismic  trace  regarded  as  a  time  series.  f'(t)  is  the  first  difference  and 
is  a  weighting  constant  dependent  on  the  sample  rate  and  noise  characteristics  at  the  station. 

The  initial  long-term  average  (LTA)  of  the  E  function  is  calculated  over  200  points  near  the 
beginning  of  the  seismogram.  It  is  then  continuously  updated  with  the  short-term  average 
(STA)  as  long  as  an  event  has  not  been  detected.  A  2-point  STA  is  used. 

The  signal  onset  is  recognized  as  the  first  time  the  ratio  STA  to  LTA  exceeds  a  threshold 
(~5).  The  algorithm  then  checks  that  the  signal  lasts  at  least  4.5  s  and  contains  at  least  4  large 
zero  crossings  before  deciding  that  a  signal  has  been  found.  If  no  signal  is  found  or  if  the  STA/ 
LTA  is  1  ‘ss  than  10,  the  trace  is  phase-free  filtered  and  re-examined.  (This  avoids  phase-free 
filtering  strong  signals.’  If  the  STA/LTA  is  less  than  5,  the  trigger  threshold  is  reduced  to  4 
and  the  trace  is  re-examined.  (This  is  to  allow  emergent  signals  to  be  picked  earlier  than  it 
would  with  a  higher  threshold.) 

After  the  trace  is  searched,  once  if  a  strong  signal  is  found  or  twice  for  a  weak  signal,  the 
program  estimates  the  position  of  the  first  arrival.  This  is  done  by  backing  along  the  trace 
from  the  point  where  the  program  triggered  to  the  first  preceding  inflection  point.  Here,  a 
marker  is  put  on  the  trace  and  the  corresponding  time  is  calculated. 

The  program's  estimate  of  the  quality  of  its  pick  is  directly  related  to  the  power  of  the  sig¬ 
nal.  If  the  STA/LTA  is  greater  than  10,  the  first  arrival  is  designated  "iP."  If  the  ratio  is  be¬ 
tween  5  and  10,  it  is  called  "eP."  A  first  arrival  with  a  ratio  less  than  5.0  is  considered  poorly 
known  and  is  labeled  "e(P)." 

After  estimating  the  first  arrival,  the  maximum  amplitude  in  the  first  5  s  is  located  and  the 
period  calculated.  The  amplitude  is  measured  from  the  raw  trace  and  converted  to  displace¬ 
ment  in  millimicrons. 

Finally,  all  the  information  obtained  is  written  to  a  summary  file,  an  example  of  which  is 
3hown  in  Table  II -7. 

Preliminary  tests  of  the  program  have  been  run  on  three  databases.  One  database  com¬ 
prises  230  seismograms,  all  the  detections  from  all  the  SRO  stations  on  3  April  1980.  From 
these  detections  the  program  picked  127  arrivals,  and  a  human  picked  141  arrivals.  The  pro¬ 
gram  rejected  7  other  seismograms  because  they  were  too  short  to  process.  Table  II— 5  shows 
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TABLE  11-5 

TEST  RESULTS  ON  3  APRIL  1980  DATABASE 

Machine  De$lgnated  Categories 

Time  Difference 

(percent) 

(Human  —  Machine) 

IP 

eP 

e(P) 

<0.5  s 

76 

48.4 

34.8 

<1.0$ 

90 

77.4 

50.0 

>2.0$ 

6 

16.1 

41.3 

Total  number  of  pick$ 

50 

31.0 

46.0 

Total  number  of  machine  pick$ 

127 

Total  number  of  machine  po$$!ble  pick$ 

134 
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the  results  of  a  comparison  between  the  machine  picks  and  the  human  picks.  The  program  made 
a  pick  on  127  of  134,  or  94.8  percent,  of  all  possible  first  arrivals.  The  human  and  machine 
picks  agreed  to  less  than  1  s  on  92  picks,  and  disagreed  by  more  than  2  s  on  27  picks.  The 
program  designated  19  of  these  27  picks  as  e(F). 

A  second  database  used  to  test  the  program  consists  of  all  the  arrivals  at  MAJO  over  a 
10-day  period  for  which  epicenters  are  available  from  the  JMSA,  PDE,  and  SDAC  bulletins. 

The  database  contains  65  seismograms  and  65  hand-picked  first  arrivals.  Results  of  a  program 
test  run  are  listed  in  Table  II-6. 


TABLE  11-6 

TEST  RESULTS  ON  MAJO  DATABASE 


- - - 

Machine  Designated  Categories 
(percent) 

(Human  —  Machine) 

iP 

eP 

e(P) 

<0.5$ 

65.0 

66.7 

33.3 

<  1 . 0  s 

92.5 

66.7 

58.3 

>2. Os 

5,0 

16.7 

16.7 

Total  number  of  picks 

40.0 

12.0 

12.0 

Total  number  of  machine  picks  64 
Total  number  of  picks  possible  65 


The  program  picked  a  first  arrival  on  98.5  percent  of  the  seismograms.  The  human  and 
machine  picks  agreed  to  less  than  1  s  on  52  of  65  possible  picks,  or  80  percent  of  the  time. 

They  disagreed  by  more  than  2  s  on  6  picks,  or  9.2  percent  of  the  time.  Two  of  these  six  picks 
including  the  greatest  disagreement  (14.5  s)  were  designated  e(P). 

The  third  database  examined,  a  set  of  25  detections  from  ANMO  on  4  October  1980,  is  part 
of  the  IDCK.  The  seismograms,  shown  in  Fig.  II - 5,  are  detection  windows  provided  by  the  al¬ 
gorithm  discussed  in  Secs.  C  and  D.  In  the  figure,  X  is  the  program's  estimate  of  the  first  ar¬ 
rival;  U  is  the  pick  provided  by  an  NEIS  analyst;  0  is  the  detection  time;  A  is  the  greatest  peak 
(or  trough)  in  the  first  5  s  of  the  signal;  and  B  is  the  adjacent  trough  (or  peak)  used  to  calculate 
ground  motion  and  period. 

Table  II— 7  is  the  output  of  the  program  for  this  set  of  seismograms.  There  is  no  informa¬ 
tion  for  seismograms  1,  5,  6,  7  or  10  since  the  program  did  not  identify  them  as  seismic  sig¬ 
nals.  The  column  of  numbers  following  the  onset  times  is  the  set  of  time  differences  between 
human-  and  machine-picked  times  (human  -  machine).  The  residual  (-79.1  s)  on  seismogram  9 
arises  from  the  windowed  seismogram  not  containing  the  event  onset  time. 

An  examination  of  the  data  reveals  that  the  machine  picks  are  probably  more  correct  than 
the  analyst  picks  on  seismograms  3,  14,  and  16  while  the  analyst  picks  are  better  on  8  and  13. 
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TABLE  11-7 

EXAMPLE  OF  PROGRAM  OUTPUT 


Seismogram 

Station 

Confidence 

Date 

Arrival  Time 

Difference  t 

Ground 

Motion 

(mp) 

Period 

(0 

Distance 

2 

ANMO 

ePD 

10/4/80 

4)21:50.1 

+0.1 

4.39 

0.95 

TELESEISM 

3 

ANMO 

e(P)D 

10/4/80 

4:44:31.2 

+0,8 

1.39 

0.65 

TELESEISM 

4 

ANMO 

iPC 

10/4/80 

4:50:11.5 

-1.2 

27.92 

1.70 

TELESEISM 

8 

ANMO 

•(P)C 

10/4/80 

14:  3:56.9 

-2.4 

1.21 

0.90 

TELESEISM 

9 

ANMO 

•<P)D 

10/4/80 

16:41:39.7 

-79.1 

3.50 

0.95 

- 

11 

ANMO 

e(P)C 

10/4/80 

18:  4:15.1 

+0.4 

0.94 

0.35 

LOCAL 

12 

ANMO 

e(P)D 

10/4/80 

18:14:41.7 

-1.9 

0.34 

0.40 

TELESEISM 

13 

ANMO 

e(P)D 

10/4/80 

19:  2:46.2 

-2.2 

3.59 

0.55 

LOCAL 

14 

ANMO 

e(P)D 

10/4/80 

19:34:49.5 

+3.7 

0.88 

0.55 

LOCAL 

15 

ANMO 

•(P)C 

10/4/80 

19:57:36.6 

- 

1.57 

0.90 

TELESEISM 

16 

ANMO 

ePD 

10/4/80 

20:16:  6.5 

+1.3 

2.44 

0.30 

LOCAL 

17 

ANMO 

•(P)D 

10/4/80 

20:28:20.8 

- 

1.05 

0.40 

LOCAL 

18 

ANMO 

*(P)D 

10/4/80 

21:  9:24.1 

- 

0.87 

0.25 

LOCAL 

19 

ANMO 

e(P)C 

10/4/80 

21:34:16.2 

- 

0.95 

0.25 

LOCAL 

20 

ANMO 

e(P)C 

10/4/80 

21:36:26.1 

- 

0.46 

0.40 

LOCAL 

21 

ANMO 

ePD 

10/4/80 

21:42:46.6 

-1.6 

1 .33 

0.60 

- 

22 

ANMO 

e(P)C 

10/4/80 

21:47:35.8 

- 

0.86 

0.45 

LOCAL 

23 

ANMO 

ePC 

10/4/80 

22:  5:  4.2 

- 

4.03 

0.50 

LOCAL 

24 

ANMO 

ePC 

10/4/80 

22:10:33.6 

- 

3.54 

0.45 

LOCAL 

25 

ANMO 

iPC 

10/4/80 

23:47:57.2 

-0.2 

11.03 

1.20 

TELESEISM 

fTime  difference  between  human  and  machine  picks  (human 

-  machine). 
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On  seismogram  4,  the  analyst  pick  is  probably  a  few  tenths  of  a  second  early,  and  the  machine 
pick  about  0.5  n  late.  It  is  difficult  to  judge  between  the  two  on  seismograms  11  and  12.  Both 
picks  may  be  several  seconds  late  on  seismogram  21.  The  events  on  17,  18,  19,  20,  and  22 
may  be  caused  by  machinery  near  the  station.  Seismograms  23  and  24  contain  eventB  timed  by 
the  program  but  missed  by  the  analyst.  The  program's  pick  appears  about  1  s  late  on  seismo¬ 
gram  23  and  about  0.5  s  late  on  seismogram  24.  The  detector  onset  time,  which  identifies  only 
the  time  that  the  detector  triggered  and  is  not  claimed  to  be  the  true  onset,  is  usually  inferior 
to  that  provided  by  this  algorithm.  Only  for  seismogram  13  is  it  better. 

In  summary,  the  program  can,  at  present,  provide  an  onset  time  within  1  s  of  what  an  ana¬ 
lyst  would  pick  about  90  percent  of  the  time  for  events  it  designates  as  iP  and  about  60  to  80  per¬ 
cent  of  the  time  for  events  it  calls  eP.  This  is  sufficient  for  the  purpose  of  event  association 
and  preliminary  location  and  should  greatly  ease  the  burden  of  the  analyst.  Furthermore,  since 
it  applies  the  same  objective  criteria  to  evaluating  the  arrivals  on  all  seismograms,  its  picks 
should  be  more  consistent  than  those  of  the  analyst. 

Research  to  improve  the  performance  of  the  program  is  proceeding  in  three  areas.  First, 
signals  so  emergent  as  to  be  impossible  for  a  human  analyst  to  time  more  accurately  than  sev¬ 
eral  seconds  sometimes  elude  the  program.  A  possible  solution  to  this  problem  is  to  alter  the 
up  ing  procedure  for  E  (long-term  average).  Second,  instead  of  merely  backing  along  the 
trace  to  the  first  inflection  point  preceding  the  trigger  point,  a  more  thoughtful  way  of  locating 
the  exact  position  on  the  first  arrival  is  being  sought.  Finally,  the  criterion  for  deciding  if  an 
event  is  local  or  teleseismic  is  being  studied  to  make  sure  of  its  applicability  at  stations  other 


than  MAJO. 


M.  W.  Shields 


F.  USE  O ,/  THE  MULTIPLE-ARRIVAL  RECOGNITION  SYSTEM  (MARS) 

FOR  TIMING  AND  IDENTIFYING  SEISMIC  BODY  WAVES 

The  most  general  waveform  containing  body  waves  also  contains  Earth  noise,  source- 
gene  if,  surface  waves  at  local  and  regional  distances,  and  randomly  scattered  body  waves 
genera,  ty  velocity  inhomogeneities  along  ray  paths  and  topographic  irregularities  at  velocity 
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disconti  ,es.  Masso  et  al.  have  described  a  method  of  waveform  processing  that  exploits 
the  nondispersive  property  of  body  waves  to  separate  them  from  the  dispersive  components  of 
a  waveform.  The  MARS  applies  a  series  of  Gaussian  filters  F  having  different  center  frequen¬ 
cies  fk  to  8  r  (nal  S,  where 


The  parameter  a  is  related  to  the  center  frequency  and  Q  of  the  filter  by 


fn2Q2(fk) 


(II —2) 


The  envelope  function  of  the  filter  output  in  the  time  domain  can  be  formed  from  the  output  spec- 

Q 

trum  and  its  Hilbert  transform  (Dziewonski  et  al.  ).  The  Q  of  the  narrowband  filters  is  chosen 
to  compromise  between  frequency  resolution  of  the  filters  and  time  resolution  of  the  envelope 
function.  Peaks  in  the  envelope  function  for  each  narrowband  filter  are  then  displayed  in  a  group 
arrival-frequency  (tg-f)  plane.  Whenever  a  sufficient  number  of  peaks  line  up  within  a  narrow 
time  band,  a  nondispersed  (body  wave)  is  detected. 


Figure  II -6  shows  the  t  -f  plane  for  a  presumed  quarry  blast  recorded  by  the  vertical  com¬ 
ponent  of  the  NSS  station  CPO.  Before  application  of  the  narrowband  filters,  the  spectrum  of 
the  waveform  is  corrected  for  the  velocity  response  of  the  short-period  band  of  an  NSS  response. 
The  waveform  is  corrected  for  velocity  rather  than  displacement  in  order  to  emphasize  the  seis¬ 
mic  energy  radiated  near  the  corner  frequency  of  body  waves.  (An  optimal  spectral  correction 
should  generally  be  based  upon  body-wave  Bpectral  characteristics  at  a  given  distance  as  well 
as  the  instrument  response.)  Q  is  taken  to  be  S.5,  giving  an  a  of  50  at  i  Hz.  The  theoreti¬ 
cal  time  uncertainty  of  the  envelope  function  at  a  center  frequency  f  is  given  by 

Thus,  at  2  Hz  At  ^  0.3  s.  The  threshold  at  which  a  peak  in  a  time  envelope  should  be  considered 
significant  should  be  based  upon  the  level  of  noise  in  the  envelope.  Here,  the  threshold  is  simply 
taken  to  be  twice  the  LTA,  The  lineup  of  envelope  peaks  in  the  t  -f  plane  shown  in  Fig.  II-6 
clearly  reveals  the  Pn,  Sn,  and  S  arrivals  as  well  as  some  unidentified  body  waves.  The  first 
three  body  waves  following  Pn  may  be  associated  with  the  P*,  Pg,  and  P  phases.  From  the  Pn, 
Sn,  and  S  phase  identifications  the  distance  is  estimated  to  be  244  *  10  km.  Note  that  the  Pn 
arrival  has  some  dispersive  character  visible  in  both  the  seismogram  and  in  the  t  -f  plane 
The  apparent  dispersion  of  Pn  may  reflect  its  character  as  an  interference  head  wave  along  the 
Moho,  a  ray  that  penetrates  below  the  Moho  forming  the  first  arrival  with  lower  relative  fre¬ 
quency  content,  and  an  infinite  series  of  multiply  reflected  and  scattered  rays  along  the  under¬ 
side  of  the  Moho  forming  the  higher-frequency  coda. 

Preliminary  results  of  applying  MARS  processing  to  teleseismic  and  regional  data  indicate 
it  to  be  a  valuable  aid  to  the  seismic  analyst  when  automatic  detections  are  refined  and  reviewed. 
The  value  of  MARS  processing  lies  in  its  ability  to  detect  arrivals  in  the  disturbed  portion  of  the 
waveform  following  a  strong  arrival.  By  aiding  the  identification  of  multiple  body-wave  arrivals, 
processing  improves  distance  estimates  of  regional  data  recorded  at  a  single  station  and  im¬ 
proves  detection  of  depth  phases  associated  with  regional  and  teleseismic  body  waves. 

V.  F.  Cormier 


G.  SYNTACTIC  ANALYSIS  OF  REGIONAL  SEISMIC  SIGNALS 

Although  many  signal-processing  techniques  are  used  in  seismic  analysis,  application  of 
these  techniques  still  requires  human  Intervention.  For  example,  to  compute  the  group -velocity 
curve  from  a  surface-wave  seismogram,  a  seismologist  must  first  identify  the  interval  contain¬ 
ing  the  surface  wave. 

In  the  detection  of  seismic  signals,  attempts  to  use  advanced  techniques  -  such  as  polariza¬ 
tion  filtering  or  prediction-error  filtering  -  have  not  been  shown  to  be  superior  to  simpler  meth¬ 
ods  because,  although  these  methods  are  powerful,  they  require  the  skilled  eye  of  a  seismologist 
to  be  applied  and  interpreted  properly.  Thus,  if  an  automatic  system  is  to  use  these  methods  to 
advantage,  it  must  use  many  of  these  same  seismological  skills. 

Most  detectors  make  one  pass  over  the  data,  declaring  a  detection  if  a  simple  threshold 
test  is  satisfied.  Unfortunately,  seismic  signals  are  often  ambiguous  and  cannot  always  be 
correctly  identified  by  a  simple  one-pass  detector.  A  better  approach  is  to  analyze  the  entire 
waveform  and  identify  its  prominent  structures  before  making  a  decision  about  it,  much  as  a 
seismologist  does. 
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As  a  first  step  toward  this,  a  computer  laboratory  known  as  SEISMQGRAMMAR  is  under 
development.  This  laboratory  is  being  used  to  investigate  techniques  for  syntactic  structural 
analysis  of  seismic  waveforms.  The  structural  analysis  of  a  waveform  is  analogous  to  the  anal* 
ysia  performed  by  a  compiler  on  a  statement  of  a  computer  programming  language.  The  wave* 
form  is  broken  up  into  features  of  morphological  significance  or  "words."  Such  features  might 
include  peaks,  inflection  points,  or  segments  of  constant  slope.  Grammar  rules  are  then  given 
that  describe  how  these  words  may  be  combined  to  form  "sentences."  Simple  descriptions  can 
be  combined  to  form  more  complicated  descriptions. 

Although  the  analogy  between  computer  language  analysis  and  waveform  analysis  is  com* 
pelllng,  the  two  domains  are  distinctly  different.  The  techniques  used  in  computer  language 
analysis  cannot  be  transferred  to  the  waveform  domain  directly.  There  are  several  reasons 
for  this,  namely: 

(1)  In  the  waveform  domain,  noise  is  always  present  and  any  parsing  scheme 
must  be  tolerant  of  it.  Thus,  one  important  issue  is  how  to  identify  the 
significant  variations  in  the  signal.  This  is  not  always  easy  since  sig¬ 
nificant  variations  in  a  signal  can  be  small,  and  the  signal  may  have 
large  variations  that  are  not  meaningful. 

(2)  Because  of  this,  the  recognition  of  words,  or  morphs,  is  often  ambiguous. 

This  can  lead  to  very  complicated  grammars,  or  combinations  of  gram¬ 
matical  and  nongrammatical  analysis. 

(3)  The  words  used  in  waveform  analysis  have  features  which  are  continuous 
variables,  such  as  slope  or  radius  of  curvature.  This  tends  to  make 
syntactic  analysis  more  complicated  by  forcing  semantic  issues  to  be 
checked  earlier  than  would  be  required  in  analyzing  a  computer  language. 

(4)  The  traditional  left-right  approach  used  in  computer  language  parsing  is 

not  robust  enough  to  work  well  in  waveform  analysis.  A  better  approach 

is  to  identify  reliable  features  first,  and  use  them  to  guide  further 
9  ”1 1 

analysis. 

!5)  Since  a  language  which  adequately  describes  a  waveform,  in  a  particular 
application,  is  usually  not  known  a  priori,  a  waveform  parsing  system 
must  allow  analyzers  to  be  developed  and  tested  easily. 

Because  of  these  problems,  SEISMOGRAMMAR  is  based  on  the  Augmented  Transition  Net¬ 
work  (ATN)  formalism.*  The  ATN  approach  provides  the  grammar  designer  with  considerable 
freedom  and  power.  Grammatical  and  nongrammatical  methods  can  be  combined  in  a  number 
of  ways,  and  grammars  can  be  built,  modified,  and  tested  easily. 

The  parsing  techniques  being  developed  are  being  applied  to  the  problem  of  recognizing  the 
S-P  interval  on  local  and  regional  seismograms.  Figure  II-7  shows  the  variation  of  signal  char¬ 
acteristics  possible  over  a  distance  range  of  6°  to  10*.  Although  the  S  arrival  is  a  prominent 
peak  in  the  signal  envelope,  the  P  arrival  is  variable.  It  can  manifest  itself  either  as  a  peak 
or  as  a  step  beginning  a  relatively  flat  coda.  Occasionally,  both  the  Pn  and  Pg  arrivals  can 
be  identified.  Automatic  detectors  can  miss  a  weak  P  arrival  completely,  thus  further  com¬ 
plicating  the  extraction  of  the  S-P  interval. 


Before  syntactic  analysis,  the  seismic  signal  is  high-pass  filtered  (to  remove  low-frequency 
background  noise),  rectified,  and  averaged  into  5-s  blocks.  The  5-s  averaging  length  removes 
most  of  the  unnecessary  detail  from  the  signal,  yet  allows  distance  to  be  estimated  to  better 
than  a  degree. 

The  resulting  envelope  is  used  as  an  input  sentence  to  the  parser.  The  words  of  the  sen¬ 
tence  are  simply  the  segments  between  consecutive  points  on  the  smoothed  envelope.  An  al¬ 
ternative  would  be  to  use  a  preprocessing  step  to  segment  the  envelope  into  a  larger  set  of  words 

9  12 

of  morphological  significance.  *  The  simpler  approach  was  chosen  since  the  appropriate 
morphs  are  currently  not  known,  and  the  ATN  formalism  allows  morphs  to  be  built  easily  by 
the  grammar  designer  as  shown  below.1  ^ 

The  ATN  method  of  parsing  was  first  developed  by  Thorne,  Bratley,  and  Dewar13  and 
14 

Bobrow  and  Fraser.  Its  current  popularity  as  a  parser  for  natural  language  is  due  to 

1 c  16  10 
Woods.  ATN  parsers  were  first  applied  to  waveform  analysis  by  Lozano-Perez.  The 
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formulation  given  here  is  similar  to  one  described  by  Charniak,  Riesbeck,  and  McDermott 
except  that  extensions  for  waveform  analysis  have  been  added.  It  is  implemented  in  FRANZ 
LISP  running  under  the  UNIX  operating  system  on  a  VAX  1 1/780. 

The  basic  idea  is  to  represent  a  grammar  as  a  transition  net.  This  is  a  set  of  nodes  (or 
states)  with  arcs  between  them.  Arcs  are  labeled  according  to  their  type,  and  an  arc  may  be 
traversed  only  if  tests  appropriate  for  its  type  are  satisfied  by  the  current  input  word.  The 
arcs  may  be  augmented  with  actions  that  do  useful  functions  such  as  building  parse  trees.  Nodes 
may  be  grouped  into  subnets,  each  subnet  being  responsible  for  parsing  one  kind  of  constituent. 
These  subnets  can  be  called  recursively  by  an  arc.  Figure  11-8  gives  an  example  of  a  simple 
ATN  that  divides  an  input  waveform  into  a  list  of  rising  or  falling  segments. 

The  arc  types  are: 

WORD  The  arc  is  traversed  if  the  current  input  word  has  the  indicated 
feature.  Before  the  arc  is  traversed,  the  current  input  pointer 
is  advanced  to  the  next  word  in  the  input.  For  example,  (word  +) 
checks  that  the  current  word  has  a  positive  slope.  A  feature  is  a 
predicate  that  the  word  must  satisfy. 

TEST  The  arc  is  traversed  if  the  arbitrary  test  is  satisfied.  For  exam¬ 
ple,  (test  t)  is  a  test  which  is  always  satisfied. 

PARSE  The  indicated  subnet  is  called  recursively  and,  if  it  succeeds,  the 
arc  is  traversed.  A  subnet  can  return  a  value  that  is  usually  a 
description  of  the  constituent  it  recognized. 

DONE  This  arc  does  a  return  from  a  subnet.  The  final  action  of  the  arc  is 
returned  as  the  value  of  the  call. 

Parsing  proceeds  as  follows.  Given  an  input  sentence  and  an  ATN' grammar,  the  arcs  em¬ 
anating  from  the  initial  state  are  inspected  in  turn.  The  first  successful  arc  is  traversed,  and 
the  arcs  at  the  new  node  are  processed  similarly.  If  a  node  is  reached  for  which  no  arc  can  be 
traversed,  the  parser  backs  up  to  the  most  recent  node  that  has  an  alternative  arc  and  tries  to 
continue  the  parse.  This  ability  to  back  up  and  pursue  another  alternative  is  valuable  in  wave¬ 
form  analysis.  This  allows  the  ATN  to  use  local  and  global  information  in  its  analysis,  which 
would  be  difficult  to  do  otherwise. 


The  ATN  grammar  shown  in  Fig.  II-8  consists  of  three  subnets.  *n»e  SEG+  net  accepts  con¬ 
secutive  segments  with  positive  slope.  Arc  actions  (not  shown)  simply  return  the  beginning  and 
end  points  of  the  resulting  segment.  The  SEG-  net  works  similarly  for  negative  segments.  The 
SEGi  net  alternatively  calls  the  SEG+  and  SEG-  subnets  and  returns  a  list  of  the  returned  seg¬ 
ments.  In  Fig,  ll-9(a),  the  effect  of  applying  this  grammar  is  shown.  The  dotted  curve  is  the 
original  envelope,  and  the  solid  one  is  the  resulting  segmentation. 

One  problem  with  this  definition  of  SEG+  and  SEG-  is  that  a  small  segment  of  opposite  sign 
will  stop  them.  It  would  be  better  to  look  ahead  slightly  to  see  if  the  segment  can  be  continued 
further  on.  This  is  what  the  SEG2  grammar  in  Fig.  11-10  does.  Its  top  level  is  the  same  as  the 
previous  one.  It  simply  returns  a  list  of  segments.  However,  the  SEG+  and  SEG-  nets  have 
been  expanded.  The  words  used  by  this  grammar  are  the  segments  produced  by  the  SEGI  gram¬ 
mar.  Four  different  word  categories  are  recognised: 

S+  A  short  segment  with  positive  slope. 

S-  A  short  segment  with  negative  slope. 

BIG+  Either  a  sharply  rising  segment,  or  a  long  segment  with  positive 
slope. 

BIG-  Either  a  sharply  falling  segment,  or  a  long  segment  with  negative 
slope. 

The  SEG+  net  is  recursive.  SEG+  will  accept  either  a  BIG+  or  a  BIG+  followed  by  an  SEG4. 

It  will  also  allow  a  segment  to  begin  with  a  short  segment  of  either  sign.  Thus,  recursion  and 
backtracking  allow  the  parser  to  look  ahead  in  the  waveform  before  deciding  on  an  interpreta¬ 
tion.  The  results  of  applying  this  grammar  are  shown  in  Fig.  II-9(b).  The  representation  is 
more  compact  than  th..t  of  the  SEGI  grammar,  though  the  general  features  of  the  envelope  are 
still  preserved.  For  a  comparison,  a  25-s  running  median  filtered  version  of  the  envelope  is 
shown  in  Fig.  II-9(c). 

The  grammars  shown  so  far  have  done  a  form  of  data  compression.  They  returned  a  simple 

list  of  the  constituents  found,  without  trying  to  determine  any  higher-level  description  of  the 

seismogram,  or  extract  any  features  such  as  the  S-P  interval.  Figure  11-11  shows  one  way  of 

building  such  a  description.  The  dotted  lines  form  a  binary  tree  structure  that  represents  the 

relationships  between  the  dominant  peaks  and  valleys  in  the  envelope.  This  approach  is  similar 

1 8 

to  that  of  Ehrich  and  Foith  except  that  the  ATN  parser  was  used  to  remove  minor  peaks,  elimi¬ 
nate  noise  and  coda,  and  build  the  peak  structure. 

Although  the  grammars  presented  here  are  simple,  they  show  the  main  features  of  the  ATN 
approach  to  structural  analysis.  Too  often,  automatic  seismic  programs,  such  as  detectors, 
develop  into  large  ad  hoc  programs  which  are  difficult  to  debug  and  test.  Typically,  a  detector 
must  perform  a  number  of  complicated  checks  and  measurements  on  the  signal  before  and  after 
a  detection.  For  example,  the  detector  might  check  if  the  detected  signal  is  of  seismic  origin 
or  not  (it  could  be  a  glitch  or  a  noise  burst),  attempt  to  refine  the  origin  time,  or  measure  other 
parameters.  The  ATN  paradigm  provides  a  convenient  framework  for  developing  such  auto¬ 
matic  processing  systems.  K.  R< 

H.  RAYLEIGH-WAVE  AZIMUTH  DETECTOR 

We  have  combined  the  work  of  Frasier, ^  Frasier  and  North,21  and  North22  in  the  develop- 

23 

ment  of  a  simple  but  effective  Rayleigh-wave  detector.  Unlike  the  more  sophisticated  Smart 
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detector,  this  one  operates  in  the  time  domain  and  on  Rayleigh  waves  only.  Unlike  P-wave  de¬ 
tectors  that  are  based  on  some  function  of  an  amplitude  ratio  taken  over  short-  and  long-time 
windows,  the  Rayleigh-wave  detector  is  tuned  to  retrograde  particle  motion  from  which  the 
azimuth  of  the  propagation  path  can  be  estimated. 

The  first  step  in  the  processing  is  to  prefilter  so  as  to  isolate  a  band  of  energy  toward  the 
long-period  end  of  the  spectrum.  This  diminishes  the  effects  of  multipathing  which  is  desirable, 
but  also  diminishes  the  dispersed  character  typical  of  shorter-period  Rayleigh  waves  and,  if 
taken  to  extremes,  will  result  in  a  signal-to-noise  problem.  This  will  occur  because  noise 
levels  increase  toward  long  periods,  beyond  a  spectral  minimum  in  the  passband  from  0.02  to 
0.04  Hz.  Typical  noise  spectra  at  SRO/ASROs  are  shown  in  Fig,  11—12  and  in  a  Texas  Instru- 
ments  Report  by  Strauss  and  Weltman.  *  Most  of  the  results  discussed  in  this  and  subsequent 
reports  come  from  data  that  have  been  bandpassed  from  0.015  Hz  (66  s)  to  0.04  Hz  (25  s).  Those 
frequencies  (periods)  define  the  3-dB  points  of  a  phaseless  3-pole  Butterworth  filter.  The  spec¬ 
trum  of  the  broadband  record  shown  in  Fig,  11-13  is  an  illustration  of  data  with  significant  signal 
strength  in  this  passband.  It  should  be  pointed  out  that,  rather  than  noise  minima,  the  important 
detection  parameter  is  signal-to-noise  ratio  which  has  not  been  studied  in  a  systematic  way. 

After  prefilter,  the  data  are  polarized  and  filtered  for  elliptical  particle  motion  in  the  man- 
20  22 

ner  proposed  by  Frasier  and  North.  Briefly  stated,  elliptical  particle  motion  is  detected  by 
first  applying  a  90®  phase  advance  to  the  horizontal  components  which  for  retrograde  elliptical 
motion  makes  the  horizontals  in  phase  with  the  vertical  component.  The  motion  in  a  time  win¬ 
dow  of  fixed  length  is  fit  in  a  least-squares  sense  to  an  ellipse  for  a  fundamental-mode  Rayleigh 
wave.  This  ellipse  will  be  elongated  and  the  correlation  coefficient  between  the  vertical  and  at 
least  one  of  the  horizontal  records  should  be  high.  The  output  of  the  detector  is  the  lom-’  axis 
of  the  ellipse  multiplied  by  the  filter  coefficient  g(t)  and  projected  onto  the  vertical  .nd  horizon¬ 
tal  axes.  The  filter  coefficient  is  either  the  correlation  coefficient  pertaining  to  9  moving  time 
window  of  fixed  length21  or  an  elliptical  filter  coefficient 


g(t)  = 


1/2 


(II-4) 
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proposed  by  Frasier.  The  time-dependent  variable  e(t)  is  the  ellipticity,  eQ  is  a  cutoff  ellip- 
ticity  (Frasier  suggested  a  value  of  0.2),  and  n  is  the  order  of  the  filter  (Frasier  suggested 
3  to  6). 

Azimuth  Az(t)  is  then  the 


arc  tan  (11-5) 

where  EW(t)  and  NS(t)  are  the  horizontal  seismograms  passed  through  the  bandpass  and  particle- 
motion  filters. 

Experience  to  date  would  suggest  that  the  ellipticity  and  correlation  filters  can  be  made  to 
work  equally  well.  Table  11-8  gives  the  detections  with  both  filters  at  stations  ANMO  and  CHTO 
for  4  October  1980.  Rayleigh-wave  trains  with  average  amplitudes  as  low  as  a  few  tens  of 
nanometers  are  detected.  The  ellipticity  filter  in  this  case  gives  fewer  detections  because  an 
eQ  of  0,2  rejects  as  not  elliptical  enough  waves  that  the  correlation  filter  passes.  The  cutoff 
value  for  the  correlation  coefficient  is  0.7,  where  unity  is  perfect  correlation.  More  work  needs 
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TABLE  11-8 

RAYLEIGH-WAVE  DETECTIONS,  4  OCTOBER  1980 


ANMO 


CHTO 


t  Half-length  of  correlation  window  32  s. 
♦  eQ  =  0. 2  and  n  is  3  (see  Eq.  (11-4)]. 


to  be  done  to  find  an  optimum  bandpass  for  the  prefilter  and  optimum  values  of  eQ  and  n  for  the 
ellipticity  filter.  T.J.  Fitch 

I.  ASSOCIATION  OF  RAYLEIGH  WAVES 

Criteria  have  been  tested  for  the  association  of  Rayleigh-wave  detections  coming  from  the 
azimuth  detector  discussed  in  the  two  preceding  SATSs  (31  March  1980  and  30  September  1980). 
The  basic  criteria  are  applied  to  paired  detections  that  are  parameterized  by  pairs  of  onset 
times,  average  amplitudes,  and  azimuths.  To  be  associated,  the  ratio  of  average  amplitude 
must  lie  in  the  range  0.1  to  10.0.  Amplitude  ratios  outside  this  range  are  unlikely  to  be  con¬ 
sistent  with  the  effects  of  anelasticity  and  spherical  spreading  on  radiation  from  a  single  source. 

Next,  the  azimuths  are  used  to  define  an  epicenter  from  which  the  path  lengths  can  be  esti¬ 
mated.  Then,  consistent  with  the  onset  and  path  lengths,  apparent  origin  times  are  tested  for 
consistency.  The  test  is  successful  if  an  origin  time  lies  within  a  time  interval  appropriate  of 
group  velocities  in  the  range  from  3.2  to  4.2  km/s.  For  mantle  waves,  the  upper  end  of  the 
group-velocity  range  would  be  more  appropriate.  If  the  amplitude  ratio  and  origin  time  tests 
are  successful,  the  paired  Rayleigh  waves  are  considered  associated.  Epicenters  from  a  num¬ 
ber  of  associated  pairs  are  then  compared  to  see  if  a  particular  source  region  is  indicated. 

Two  sets  of  Rayleigh-wave  detections  have  been  analyzed  in  this  way  with  the  following  re¬ 
sets.  Table  II-9  shows  epicenters  from  azimuth  pairs  pertaining  to  an  Mg  6.1  shallow  earth¬ 
quake  in  the  Fiji  Islands  region.  Azimuths  at  CHTO,  BOCO,  ANTO,  and  BCAO  differ  by  more 
than  10°  from  the  appropriate  great-circle  azimuth.  In  particular,  the  paths  tc  CHTO  and  ANTO 
which  are  largely  tectonic  in  nature  are  expected  to  produce  multipathing  which,  in  turn,  will 
result  in  large  discrepancies  between  apparent  and  great-circle  azimuths.  The  apparent  azi¬ 
muth  at  CHTO  is  90.2°,  whereas  the  great-circle  azimuth  is  107.6°.  At  ANTO,  the  apparent 


TABLE  11-9 

ASSOCIATED  RAYLEIGH  WAVES  FROM  FIJI  ISLANDS  EARTHQUAKE 
(4  October  1980,  4:37:34.0) 

Station  Pair 

Epicenter 

Origin  Timet 
Difference 

r“ 

Latitude 

(deg) 

Longitude 

(deg) 

0)  (2) 

ANMO -GUMO 

-17.6 

178.2 

-80 

ANMO-TATO 

-21.4 

172.7 

-362 

CHTO  -GUMO 

10.6 

148.1 

-30 

GUMO  -  NWAO 

-18.6 

179.4 

-256 

NWAO-TATO 

-21.6 

173.0 

304 

ANTO  -BCAO 

-7.5 

-123.2 

-368 

PDE  epicenter 

-15.76  179.31 

fO.T.  (station  1)  — O.T.  (station  2)  where  O.T.  =  onset  time  -distance/4.0. 


TABLE  1 1-10 


ASSOCIATED  RAYLEIGH  WAVES  FROM  PRESUMED  EARTHQUAKE 
IN  WEST  AFRICAt 


Station  Pair 

Epicenter 

Origin  Time* 
Difference 

Latitude 

(d«g) 

Longitude 

(deg) 

0)  (2) 

ANMO-ANTO 

15.2 

0.1 

42 

CHTO  —  NWAO 

6.7 

8.7 

-204 

ANMO  -  BOCO 

6.9 

-3.7 

-542 

CHTO  -ANTO 

12.8 

29.0 

464 

CHTO  -BOCO 

8.9 

15.4 

366 

CHTO  —  BCAO 

7.6 

12.2 

-158 

NWAO  -  TATO 

2.1 

-141.3 

408 

NWAO  -  BCAO 

16.1 

-7.2 

18 

TATO  -BOCO 

9.0 

-156.2 

454 

TATO  -BCAO 

6.7 

34.2 

318 

ANTO  -  BCAO 

15.5 

15.8 

-54 

BOCO  -  BCAO 

9.3 

8.4 

-120 

ANTO  -KONO 

19.8 

4.6 

104 

ANTO  -KONO 

12.1 

4.1 

44 

fNo  PDE  location. 

t  O.T.  (stavion  l)-O.T.  (station  2)  where  O.T.  =  onset  time  —  distance/4,0. 


azimuth  is  36.3°  and  the  great-circle  azimuth  is  58.4°.  In  contrast,  at  ANMO  these  azimuths 
differ  by  less  than  10°.  The  difference  in  apparent  origin  times  for  each  pair  of  detections 
ranges  from  a  fraction  of  a  minute  to  more  than  6  min.  If  a  nominal  group  velocity  <f  4,0  km/s 
is  assumed,  the  origin  time  differences  correspond  to  an  uncertainty  in  epicenter  of  as  much 
as  20°.  Clearly,  the  criteria  by  which  detection  pairs  are  judged  to  be  associated  or  not  are 
inadequate  for  accurate  epicenter  determinations. 

Table  11-10  contains  epicenter  pairs  for  a  presumed  earthquake  at  a  latitude  of  about  15° 
with  a  longitude  that  is  uncertain  by  as  much  as  20'  in  arc  distance.  Most  of  the  epicenters 
lie  in  West  Africa,  which  is  consistent  with  the  largest  amplitude  recorded  at  BCAO.  However, 
a  longitude  passing  through  West  Africa  would  give  an  epicentral  distance  to  BCAO  with  the 
range  10°  to  20°  which  is  too  close  to  be  consistent  with  the  dispersed  wave  trains  shown  for 
event  2  in  Fig.  11-14.  The  dispersion  and  azimuth  are  consistent  with  an  earthquake  on  the  mid- 
Atlantic  ridge.  Although  this  association  of  Rayleigh  waves  from  this  azimuth  detector  has  not 
yet  given  much  encouragement  of  accurately  determining  epicenters  by  this  means,  a  rough 
origin  time  and  location  could  be  used  to  search  the  short-period  records  for  P-waves  that  have 
gone  undetected  in  a  first  pass  through  the  data. 

T.  J.  Fitch 


J.  SRO/ASRO  RAYLEIGH-WAVE  DETECTIONS 

The  long-period  seismograms  from  the  SRO/ASROs  contain  an  abundance  of  surface  waves 
most  of  which  cannot  be  associated  with  events  located  by  the  standard  technique  using  short- 
period  P  times.  Conversely,  only  a  fraction  of  the  events  listed  in  the  PDE  or  ISC  bulletins 
can  be  correlated  with  a  Rayleigh-wave  detection.  The  detections  considered  here  are  from 
the  azimuth  detector  previously  described  in  this  report. 

Examples  of  these  raw  and  processed  data  from  stations  ANMO  and  CHTO,  shown  in 
Figs.  II-l  5(a)  and  (b),  illustrate  some  of  the  major  problems  encountered  when  evaluating  in¬ 
dividual  detections.  A  signal  was  considered  detected  if  elliptical  particle  motion  was  cor¬ 
related  on  the  vertical  and  at  least  one  of  the  horizontal  components  to  better  than  0.7  for  more 
than  60  s.  The  filtered  records  for  ANMO  show  that  the  detector  often  switches  off  and  on  some¬ 
times  sever?,  times  during  the  passage  of  what  looks  like  a  single  wave  train  to  the  observer. 
This  problem  is  exacerbated  by  multipathing  that  is  seen  as  a  beating  pattern  on  Rayleigh  waves 
that  have  traveled  teleseismic  distances. 

Parameters  of  mean  azimuth  in  the  detection  window,  mean  amplitude,  onset  time,  and 
duration  were  stored  for  each  detection.  Table  11-11  summarizes  the  results  of  a  comparison 
between  Rayleigh-wave  detections  and  PDE  locations  for  4  October  1980.  At  most,  only  half 
of  the  located  events  have  detectable  Rayleigh  waves  at  a  given  station,  and  at  a  sensitive  site 
such  as  the  one  for  BCAO  in  Central  Africa  there  will  be  many  times  as  many  detections  as 
there  are  local  events. 

These  disparities  are  not  surprising  because  the  body  phases,  of  which  the  P-waves  are 
prime  examples,  and  the  surface  waves  have  fundamentally  different  excitation  functions  that 
are  strongly  depth  depend  jnt  in  the  case  of  fundamental -mode  surface  waves.  P-waves  are 
excited  preferentially  by  underground  explosions  or  earthquakes  at  subcrustal  depths.  Of  the 
located  events  on  4  October,  two  (m^  4.7  and  4.8)  were  placed  in  the  mantle  and,  as  expected, 
Rayleigh  waves  from  them  were  not  detected  at  teleseismic  distances.  In  contrast,  Rayleigh 
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TABLE  11-11 

PDE  LOCATIONS  FOR  4  OCTOBER  1980 

VS  RAYLEIGH-WAVE  DETECTIONS 

Sta'  ion 

Number  of 
Rayleigh-Wave 
Detections 

Number  Correlated 

With  10  Events 
of  m^  >  3.5t 

AN  MO 

11 

3 

CHTO 

15 

4 

NWAO 

13 

3 

ANTO 

20 

4 

BOCO 

8 

2 

6CAO 

62 

3 

fStation  to  epicenter  azimuth  agree  within  ±30°,  and 

Rayleigh-wave  arrival  times  predicted  from  locations 
agree  with  detector  onset  times  within  a  few  minutes, 
assuming  a  group  velocity  of  3,  8  km/s. 

waves  from  a  small  earthquake  (m^  5.2,  O.T.  03:24:49.0)  on  4  October  were  located  in  the 
Philippine  Islands  region  but  yielded  no  short-period  P-wave  detections  at  even  the  closer  of 
the  SRO/ASROs  -  CHTO,  GUMO,  TATO,  and  NWAO  -  whereas  Rayleigh  waves  with  reasonable 
onset  times  and  azimuth  were  detected  at  these  stations.  The  Rayleigh-wave  detections  can  be 
formally  associated  with  a  particular  event  if  amplitudes,  onset  times,  and  azimuths  meet  cer- 
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tain  criteria  that  are  discussed  in  Sec.  III. 


T.  J.  Fitch 
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etecti^.is  from  ANMO  SPZ  4  October  1980.  X  is  program 
ick;  U  is  analyst  pick  provided  with  IDCE  parameter  data- 
tection  time;  A  is  largest  peak  in  first  5  s. 


40 


Fig.  II -6.  The  tg-f  plane  and  seismogram  for  an  event  recorded  by  short-period 
band  of  CPO.  Start  time  is  4  October  1980,  17:07:10.45  GMT. 
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Fig.  H-7.  Variation  of  seismic  envelope  $ 
shape  between  6*  and  9*.  Original  seis-  a 
mograms  were  recorded  at  ANMO,  " 
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Fig.  11-8.  A  simple  ATN  grammar.  Subnet 
for  SEG-  (not  shown)  is  similar  to  that  of 
SEG+. 
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Fig.  II— 9.  Data  representation  using  (a)  SEG1  grammar,  and  (b)  SEG2 
grammar.  Effect  of  a  25-s  running  median  smoother  is  shown  for 
comparison  in  (c). 


Fig.  11-10.  SEG2  grammar.  Subnet  for  SEG-  (not  shown)  is  similar 
to  that  of  SEG+. 


Fig.  II— 1 1.  Hierarchical  structure  (dotted  lines)  imposed  by  relationships 
between  peaks  and  valleys  in  seismic  envelope  (solid  curve). 
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Fig.  11-12.  SRO/ASRO  long-period  noise  spectra  from  1  April  1980. 
Instrument  response  and  spectra  smoothed  by  averaging  over  seven 
to  eight  51 2 -point  time  samples.  Sampling  interval  is  4  s. 
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Fig.  11-13.  Spectra  of  a  broadoand  SRO  recording  at  ANMO.  Spectra  smoothed 
by  averaging  eight  4096 -point  time  samples.  Sampling  interval  is  i/20  s.  Notice 
spectral  lines  in  band  from  1  to  Nyquist  frequency  of  10  Hz. 
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III.  GENERAL  SEISMOLOGY 

A.  FURTHER  SURFACE-WAVE  ANALYSIS  OF  EASTERN  KAZAKH 
PRESUMED  EXPLOSIONS 

1  2 

In  the  two  preceding  SATSs,  '  we  studied  several  presumed  explosions  in  the  eastern  part 
of  the  Soviet  test  site  in  Eastern  Kazakh.  For  one  of  these  (on  7/7 /?9),  the  Rayleigh  waves  were 
reversed  in  phase  compared  with  the  others  at  all  observing  stations,  with  a  pronounced  non¬ 
circular  radiation  pattern  and  enhanced  Love-wave  radiation  pattern.  We  have  since  extended 
the  number  of  events  studied  to  those  given  in  Table  III— 1.  A  detailed  analysis  of  these  events 
(carried  out  as  before)  reveals  that,  as  well  as  the  event  on  7/7/79,  another  event  on  8/18/79 
generated  Rayleigh  waves  which  were  reversed  in  phase  at  all  stations  compared  with  the 
"normal"  events  (hereafter  called  type  I),  with  a  similar  Rayleigh-wave  radiation  pattern  and 
enhanced  Love-wave  amplitudes  (type  III).  In  addition,  a  third  event  (on  11/4/78)  is  reversed 
in  phase  at  two  stations  only  (KAAO  and  MAJO)  which  correspond  to  observation  points  in  oppos¬ 
ing  quadrants  of  the  radiation  pattern  (type  II). 


TABLE 

lll-l 

REVISED  LOCATIONS  AND  MAGNITUDES  OF  EASTERN  KAZAKH  EVENTS 

Event 

Date 

Origin  Time 
(h:m:s) 

Latitude 

E(Lat) 

(km) 

Longitude 

E(Lon) 

(km) 

mb 

M 

s 

SD(Ms) 

N(Ms) 

1 

9/15/78 

02:36:57.3 

49.920 

2.95 

78.843 

3.62 

6.0 

3.99 

0.14 

5 

2 

1 1/4/78 

05:05:57.5 

50.050 

4.18 

78.918 

3.07 

5.6 

3.75 

0.12 

7 

3 

4/23/79 

02:56:57.0 

49.896 

3.37 

78.807 

4.63 

6.3 

4.14 

0.08 

6 

4 

7/7/ 79 

03:46:57.4 

50.013 

4.58 

78.947 

4.11 

5.8 

3.95 

0.49 

5 

5 

8/4/79 

03:56:57.2 

49.887 

3.51 

78.838 

4.23 

6.1 

4.23 

0.07 

6 

6 

4/18/79 

02:51:57.3 

49.953 

3.40 

78.912 

3.34 

6.1 

3.84 

0.17 

7 

7 

IQ/24/79 

03:16:56.9 

49.987 

4.28 

78.949 

3.68 

6.0 

4.18 

0.15 

5 

8 

]  2/2/79 

04:36:57 .5 

49.903 

2.81 

78.817 

4.26 

6.0 

4.25 

0.05 

7 

9 

1 2/23/79 

04:56:57.6 

49.930 

3.05 

78.716 

3.05 

6.1 

3.96 

0.02 

4 

Date,  origin  time,  and  body-wave  magnitude  mb  from  NEIS  bulletin.  Locations  revised  as 
described  in  Ref.l,  E(Lat)  and  E(Lon)  are  standard  errors  in  latitude  and  longitude.  Surface-wave 
magnitude  M,  is  mean  of  N(M,)  observations  at  stations  of  SRO/ASRO  networks,  with  standard 
deviation  of  SD(M$) . 

As  we  have  pointed  out  previously,^  the  surface-w.ive  radiation  for  the  types  II  and  III  events 
can  be  explained  by  the  addition  of  induced  thrust  faulting  to  the  explosion  monopole.  Induced 
faulting,  though  in  this  case  of  strike-slip  type,  has  been  suggested  to  explain  observations  simi¬ 
lar  to  those  of  the  type  II  event  for  some  explosions  at  the  Nevada  Test  Site  (NTS).  Thrust  fault¬ 
ing  has  been  proposed  independently^'^  as  a  possible  mechanism  for  certain  features  of  the 
elastic -wave  radiation  from  some  underground  explosions. 
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1  2 

In  our  previous  studies  ’  we  noted  that  the  imaginary  part  of  the  far-field  complex  spectrum 
of  both  Rayleigh  and  Love  waves  is  very  poorly  excited  by  shallow  sources.  In  our  earlier  at¬ 
tempts  to  invert  surface -wave  data  from  Eastern  Kazakh,  we  carried  out  inversion  in  the  time 

5 

domain  for  the  resolvable  parts  of  the  seismic  moment  tensor  in  the  formulation  of  Mendiguren. 
We  have  subsequently  abandoned  the  time -domain  inversion  procedure  because  of  its  extreme 
sensitivity  to  small  phase  changes,  and  now  invert  the  data  in  the  frequency  domain  assuming  in¬ 
duced  thrust  faulting  and  the  explosion  monopole  to  be  the  only  contributors  to  the  surface -wave 
radiation.  This  extremely  restricted  model  ignores  some  proposed  features  of  the  seismic 
source  for  explosions  such  as  spall,  the  contribution  from  which  is  however  inextricably  linked 
with  that  of  the  monopole  in  the  surface -wave  radiation.  There  is  no  clear  indication  of  spall  in 
the  short -period  P-waves  for  any  of  the  events  studied  here,  though  this  need  not  preclude  its 
occurrence  at  very  short  intervals  after  the  explosion  itself. 

In  the  absence  of  the  imaginary  part,  the  far-field  Rayleigh-wave  radiation  for  a  shallow 
source  is  given  by 

Fr(6  )  =  ApSR(R1  -  R2  sin  29  +  R3  cos  26  ) 

.  2 

where  FR(6),  AR,  and  SR  are  as  defined  previously  (depending  additionally  upon  frequency  and 
source  depth),  6  being  the  azimuth  of  observation.  Assuming  the  source  to  consist  only  of  a 
monopole  of  seismic  moment  E  and  thrust  faulting  of  moment  Q, 

ri=4e-|q 

R2  =  Q  sin  20/2 
R3  =  Q  cos  20/2 

where  the  strike  of  the  thrust  fault  is  given  by  0.  Thus  E,  Q,  and  <p  are  given  by 
Q  =  2(R2  +  R3 )1/2 
E  =  1.5  R1  +  1.25  Q 

tan'1  (R2/R3) 

*  =  - 2 -  • 

The  data  we  wish  to  invert  are  the  real  part  of  the  source  spectrum:  we  have  observation 
only  of  the  amplitude  variations  with  azimuth  and  phase  differences  between  events.  To  obtain 
the  real  part  from  this  information,  we  make  the  following  assumptions: 

(1)  The  imaginary  part  of  the  spectrum  cannot  be  appreciably  excited  for 
shallow  sources. 

(2)  The  six  "normal"  events  (type  I)  with  the  smallest  amplitude  variation 
with  azimuth  are  such  that  the  monopole  term  dominates,  and  the  source 
initial  phase  is  it.  Observations  of  phase  reversals  at  a  given  station 
compared  with  these  normal  events  indicate  a  source  initial  phase  of  0. 

(3)  The  real  part  of  the  spectrum  is  consequently  given  by  (-amplitude  spec¬ 
trum)  for  unreversed  signals  and  by  (+amplitude  spectrum)  for  reversed 
signals. 
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We  are  forced  to  make  assumption  (2)  and  hence  (3)  by  the  lack  of  any  calibration  of  propaga¬ 
tion  effects  for  the  paths  used.  The  above  assumption  conveniently  sidesteps  the  problem  of  cal¬ 
ibrating  the  dispersion  for  each  path,  and  we  are  left  with  the  question  of  attenuation.  We  have 

-4 

chosen  an  attenuation  coefficient  y  of  1  X  10  /km  over  the  period  range  20  to  45  s  to  correct  the 

observed  amplitudes  back  to  the  source.  This  rather  arbitrary  choice  falls  between  the  largest 

6  7 

and  smallest  Q  estimates  of  Burton  and  is  toward  the  lower  end  of  the  measurements  of.  Patton. 

We  have  found  that  different  choices  of  average  y  (e.g.,  0.5  or  2.0  X  ltT^/km)  result  in  only  mar¬ 
ginal  changes  in  the  results  of  the  inversion.  For  such  disparate  paths  as  those  considered  here, 
any  gross  average  must  for  some  paths  be  substantially  in  error,  but  we  are,  again  in  the  ab¬ 
sence  of  a  calibration  event  of  known  mechanism,  unable  to  assign  different  attenuation  coeffi¬ 
cients  to  each  path.  A  very  large  difference  exists  between  both  spectral  and  time-domain  am¬ 
plitudes  at  stations  CHTO  and  SHIO,  though  these  differ  only  slightly  in  both  azimuth  and  distance 
from  the  source.  We  have  somewhat  arbitrarily  decided  to  multiply  amplitudes  at  CHTO  by  3  in 
order  to  bring  them  into  closer  agreement  with  those  at  SHIO.  We  believe  the  amplitudes  at 
CHTO  to  be  low  due  to  defocusing  of  the  ray  paths  by  oblique  incidence  at  the  high  velocity- 
contrast  boundary  between  the  Tibetan  plateau  and  the  structure  to  its  south,  and  ray-tracing 

g 

experiments  in  laterally  varying  velocity  models  of  this  area  demonstrate  such  an  effect.  It 
should  be  noted  here  that  the  effects  of  geometrical  spreading  over  the  surface  of  the  globe  are 

incorporated  in  the  term  SD  given  above. 

K  5  9 

In  the  moment  tensor  inversion  schemes  as  implemented  by  Mendiguren  and  Patton,  the 

data  to  be  inverted  consist  of  the  real  and  imaginary  parts  of  the  spectrum  at  a  range  of  frequen¬ 
cies  as  observed  over  a  number  of  stations  or  azimuths.  In  general,  this  is  necessary  because 
of  the  differing  variation  of  the  excitation  functions  AR,  B^,  and  CR  with  frequency.  This  is  not 
strictly  necessary  for  very  shallow  sources  since  the  spectra  are  determined  in  this  case  only 
by  Ap,  scaled  by  a  factor  (which  may  be  negative)  depending  on  the  value  of  Mij  and  the  azimuth. 
Thus,  provided  the  propagation  corrections  and  various  other  basic  assumptions  such  as  that  of 
a  step  source  time  function  are  correct,  the  same  solution  would  be  provided  by  observations  at 
a  single  frequency  only.  We  may  thus  choose  to  invert  the  data  in  several  different  forms: 

(1)  At  a  single  frequency  point  where  we  believe  that  variations  in  propagation 
effects  such  as  spreading  and  attenuation  are  smallest,  such  as  the  longest 
periods  available  in  the  spectrum. 

(2)  At  a  number  of  frequencies,  taking  advantage  of  the  redundancy  thus  pro¬ 
vided  to  effectively  average  several  estimates. 

(3)  Integrate  the  spectrum  over  a  certain  frequency  range,  the  excitation  func¬ 
tion  then  being  the  corresponding  integral  of  A^S^. 

After  several  experiments,  we  settled  upon  method  (3)  above,  integrating  over  the  period 
range  20  to  45  s,  as  providing  the  most  stable  data  for  inversion,  probably  because  it  averaged 
out  variations  which  could  not  otherwise  be  removed  without  excessive  smoothing.  Figure  I II - 1 
shows  this  form  of  the  input  data  for  the  inversion,  consisting  of  integrals  of  the  real  part  over 
the  chosen  frequency  range.  Negative  values  indicate  initial  phase  of  jr  (i.e.,  for  "normal"  phase 
seismograms),  and  positive  values  indicate  an  initial  source  phase  of  0  (reversed  compared  with 
normal  seismograms). 


The  results  of  a  least-squares  inversion  of  the  data  shown  in  Fig.  Ill- 1  are  given  in 
Table  III-2,  where  R^,  R2>  and  R3  are  the  e-invariant,  sin  (26),  and  cos  (26)  coefficients  pro¬ 
viding  the  best  fit  to  the  data.  The  number  of  stations  used  in  each  inversion  and  the  ratio  of 
the  residual  rms  to  that  of  the  data, itself  are  also  given.  The  latter  is  smallest  for  the  type  I 
events  exhibiting  the  least  azimuthal  variation,  and  larger  for  the  types  II  and  III  events.  We 
discovered  that  the  large  residual  rms  for  event  2  could  be  considerably  reduced  (by  50  percent) 
by  changing  the  sign  of  the  real  part  (i.e.,  reversing  the  phase)  at  station  ANTO,  but  such  a  step 
is  not  warranted  by  the  phase  difference  data.  We  have,  however,  no  real  explanation  for  the 
greater  difficulty  in  obtaining  a  fit,  as  measured  in  terms  of  residual  rms,  for  the  types  II  and 
III  events.  Misfits  were  largest  for  the  stations  near  the  "nodes,"  i.e.,  of  smallest  amplitude: 
this  is  a  common  feature  of  inversions  of  radiation  patterns  since  amplitudes  in  these  directions 
are  generally  overestimated. 

From  the  results  of  the  inversion  we  have  deduced  the  seismic  moments  of  the  explosion 
monopole  and  double -couple  components,  assuming  the  latter  induced  faulting  to  be  of  thrust  type, 
and  these,  together  with  the  ratio  of  the  moment  of  faulting  to  explosion  (or  F-factor),  are  given 
in  Table  III-2.  As  we  could  have  deduced  from  Fig.  Ill -1 ,  F  is  smallest  for  the  type  I  events  and 
largest  for  those  of  type  III,  for  which  F  exceeds  2.0,  as  required  for  phase  reversal  at  all  azi¬ 
muths.  We  do  not  really  understand  the  very  large  F  (13.28)  obtained  for  event  6,  which  is  pri¬ 
marily  due  to  the  exceptionally  low  value  obtained  for  the  explosive  moment.  Excluding  this,  we 

22 

see  that  the  values  of  the  moment  of  the  explosion  range  only  from  9.06  to  20.38  x  10  dyn-cm, 

reflecting  the  small  range  of  explosion  yields  which  can  be  deduced  from  the  similarly  small 

10 

range  of  mfa  values  for  these  events.  Blandford  has  obtained  an  expression  for  corner 


RESULTS  OF 

TABLE  111—2 

RAYLEIGH-WAVE  INVERSION 

Event 

R1 

R2 

R3 

N 

Fit 

Q 

E 

F 

Azimuth 

1 

5.30 

0.84 

-1.30 

6 

0.21 

3.09 

11,81 

0.26 

73.5 

2 

0.76 

3.07 

0.80 

8 

0.49 

6.34 

9.06 

0.79 

37.7 

3 

6.39 

2.76 

1.96 

8 

0.15 

5.60 

16.59 

0.34 

49,9 

4 

-10.60 

9.84 

-0.48 

8 

0.30 

21.16 

10.55 

2.01 

34.2 

5 

7.68 

2.99 

3.88 

8 

0.17 

7.09 

20.38 

0.35 

61.2 

6 

-5.16 

2.83 

-1.90 

9 

0.28 

6.59 

0.49 

13.28 

29.6 

7 

7.16 

2.75 

1.69 

7 

0.04 

6.50 

18.86 

0.34 

28.9 

8 

9.58 

0.32 

1.73 

9 

0.21 

4.43 

19.91 

0.22 

85.8 

9 

4.23 

0.65 

-1.59 

8 

0.07 

3.42 

10.63 

0.32 

78.9 

R],  R2,  and  R3  are  radiation  pattern  coefficients  as  defined  in  text  (units  10^  dyn-cm),  N  is  the 
number  of  stations  used,  and  the  fit  is  measured  as  the  ratio  of  the  residual  rms  to  the  data  rms.  Assum¬ 
ing  induced  thrust  faulting,  the  respective  moments  of  the  thrust  fault  and  explosion  monopole  are  given 
by  Q  and  E  (units  10^  dyn-cm)  and  their  ratio  F  =  Q/E  is  listed  in  the  next  column.  The  final 
column  gives  the  strike  (measured  in  degrees  West  of  North)  of  the  assumed  thrust  fault  contribution. 
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frequency  as  a  function  of  yield  from  data  for  NTS  explosions;  for  100-kT  yield  his  formula 
gives  a  corner  frequency  of  0.8  Hz  which,  assuming  omega-square  scaling,  implies  a  DC  seis- 
mic  moment  of  ~10  dyn-cm,  in  good  agreement  with  the  values  we  have  obtained.  It  is  worth 
emphasizing  that  a  simple  model,  such  as  an  omega-square  model  with  such  a  DC  value  of  mo¬ 
ment,  is  incapable  of  generating  the  large  short-period  P-waves  generated  by  these  explosions; 
and  a  more  complicated  source  model,  such  as  that  of  Aki  et  al. } 1  is  required  to  explain  the 
large  differences  in  the  excitation  of  high-  and  low-frequency  for  explosions  which  leads  to  the 

M„:m.  discriminant, 
s  b 

In  Fig.  III-2(a-b)  we  show  two  comparisons  of  the  monopole  explosion  moment  with  short- 
period  body-wave  amplitudes.  Figure  III-2(a)  compares  the  average  amplitude,  obtained  from 
log(mb-5)  where  mfa  is  that  given  in  the  PDE,  with  the  explosion  moment.  This  is  somewhat 
unconvincing,  primarily,  we  believe,  because  of  the  large  variations  in  amplitude  from  station 
to  station  as  exemplified  in  the  short-period  P-waves  observed  at  the  SRO/ASRO  sites.  To  at¬ 
tempt  to  obtain  a  more  consistent  measure  of  short-period  amplitudes  we  have,  for  each  of  the 
SRO/ASRO  stations,  calculated  the  ratio  of  the  P-wave  maximum  amplitude  to  that  for  event  8, 
and  averaged  this  for  the  same  set  of  stations  for  each  event.  The  results  are  shown  in 
Fig.  Ill -2 (b)  where  a  slightly  more  convincing  correlation  is  obtained.  By  both  measures,  how¬ 
ever,  the  explosive  moment  obtained  for  event  6  is  excessively  low. 

We  have  investigated  the  solution  for  event  6  in  some  detail  and  found  no  errors  in  our  re¬ 
duction  of  the  data  for  this  event.  The  data  from  KONO  are  a  crucial  factor  in  the  solution.  By 
reducing  the  spectral  amplitude  given  in  Fig.  III-l  for  KONO  by  two-thirds,  an  apparently  more 
reasonable  solution  with  F  ~  3  can  be  obtained,  but  there  is  no  reason  to  assume  that  the  data 
for  this  station  are  in  error  and,  in  fact,  the  spectral  amplitude  variations  shown  for  KONO  for 
each  event  in  Fig.  Ill- 1  exactly  mirror  the  time-domain  measurements.  The  spall  or  slapdown 
phase  may  conceivably  have  been  unusually  large  for  this  event  and  largely  canceled  out  the  mono¬ 
pole  contribution  in  the  long-period  part  of  the  spectrum.  No  indication  can  be  seen  of  such  a 
phase  in  the  short -period  data  for  this  event,  but  this  does  not  necessarily  mean  that  spall  is  not 
a  significant  factor.  It  is  unfortunately  inseparable  from  the  monopole  contribution  in  the  long- 
period  surface-wave  spectra.  We  have  assumed  step-function  time  dependence  for  both  monopole 
and  faulting  contributions.  That  this  is  not  strictly  valid  is  shown  by  the  well-known  observation 
that  the  dominant  period  of  the  maximum  amplitude  surface  waves,  from  larger  nuclear  explosions 

at  least,  is  much  shorter  (8  to  12  s)  than  that  for  earthquakes  at  similar  distances  as  seen  on  the 
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WWSSN  instrument,  though  not  on  the  narrower-band  SRO/ASRO  systems.  Several  studies 
have  indicated  a  different  source  time  function  for  explosions,  which  is  apparently  yield-dependent. 

The  deduced  seismic  moment  of  the  inferred  thrust  faulting  shows  much  more  variation  than 
that  for  the  explosive  component.  A  possible  check  on  these  values  is  provided  by  the  Love  waves, 
which  are  of  course  generated  by  the  faulting  component  only.  We  have  integrated  the  Love- 
wave  spectra  of  each  seismogram  over  the  period  range  20  to  45  s,  as  we  did  previously  for  the 
Rayleigh-wave  data,  correcting  for  spreading  and  an  attenuation  of  2  x  10  ^/km  (twice  that  used 
for  the  Rayleigh  waves).  An  average  value  of  this  was  determined  for  each  event,  and  the  results 
are  plotted  in  Fig.  III-3  against  the  predicted  maximum  value  of  this  amplitude  integral,  given 
simply  by 

HALSL  • 


51 


We  aee  that  the  amplitudes  of  the  observed  Love  waves  are  well  predicted  by  the  values 
obtained  for  the  moment  Q  of  the  (thrust)  double-couple  component.  We  have  not  attempted  to 
compare  observed  and  theoretical  amplitude  radiation  patterns  since  we  do  not  believe  we  have 
sufficient  azimuthal  resolution  to  do  so. 

The  strike  of  the  assumed  induced  thrust  faulting,  given  in  Table  III-2,  is  fairly  restricted, 

ranging  from  29°  to  80°  West  of  North.  If  we  restrict  ourselves  to  events  with  high  F-factors 

(events  2,  4,  and  6),  for  which  the  faulting  contribution  is  strongest,  the  strike  varies  only  from 

29°  to  38°  West  of  North.  While  we  are  not  aware  of  the  detailed  tectonics  of  this  region  of 

Eastern  Kazakh,  it  is  interesting  to  note  that  this  is  very  close  to  the  strike  (40°  W  of  N)  of  a 

14 

large  fault  clearly  visible  on  Landsat  imagery  of  this  general  region,  located  however  some 
70  km  to  the  west  and  additionally  apparently  of  strike-slip  nature.  The  consistency  of  the  fault 
azimuths  obtained  indicates  some  fairly  large-scale  regional  tectonic  influence  governing  the 

3 

orientation' of  the  double -couple  component  of  the  source.  In  this  respect  we  differ  from  Masse 
who,  while  advocating  thrust  faulting  as  a  significant  contribution  to  seismic  radiation  from  ex¬ 
plosions,  denies  that  this  represents  "induced"  faulting  governed  by  tectonics. 

We  believe  that,  within  the  limitations  of  the  source  model  we  have  assumed  (point  source, 
step-function  time  dependence  for  both  explosive  and  faulting  contribution,  and  no  consideration 
of  spall),  we  have  been  fairly  successful  in  explaining  most  features  of  the  surface-wave  radiation 
from  the  9  explosions  we  have  studied  here.  That  this  model  is  not  entirely  adequate  is  illus¬ 
trated  by  the  ludicrously  high  F-factor  obtained  for  event  6,  but,  nevertheless,  we  feel  that  we 
have  demonstrated  that  induced  thrust  faulting  is  a  major  contributor  to  the  surface -wave  radia¬ 
tion  from  Eastern  Kazakh  explosions.  ^  G  „ 

t!  J.'  Fitch 


B.  INDUCED  THRUST  FAULTING  AT  NTS? 

In  Sec.  A  above,  we  indicated  that  Rayleigh-wave  observations  of  Eastern  Kazakh  explosions 
in  a  small  area  to  the  north  of  the  Shagan  River  can  be  adequately  explained  by  a  source  model 
combining  an  explosion  monopole  and  a  double -couple  component  of  thrust  faulting  type.  The 
amount  of  such  induced  faulting  can,  when  sufficiently  large,  cause  phase  reversals  at  some  or 
all  azimuths.  The  moments  obtained  for  the  monopole  and  thrust  faulting  components  of  the  com¬ 
bined  source  are  well  correlated  with  short-period  P-wave  amplitudes  and  long-period  Love- 
wave  amplitudes,  respectively.  These  conclusions  are  in  good  agreement  with  an  explosion 
source  model  suggested  by  Masse.3 

Observations  of  noncircular  radiation  patterns  and  phase  reversals  in  some  quadrants  for 

some  NTS  explosions  may  also  be  explained  by  such  induced  thrust  faulting,  rather  than  the 

14 

strike-slip  faulting  which  was  assumed  by  Toksoz  and  Kehrer,  with  little  justification  from  other 
evidence  such  as  tectonic  strain  release  patterns.  We  prefer,  on  purely  physical  grounds,  the  in¬ 
ducement  of  thrust  or  normal  faulting  over  that  of  strike-slip  faulting,  since  it  seems  to  us  more 
reasonable  that  vertical,  rather  than  horizontal,  displacements  are  more  likely  to  occur  as  a  re¬ 
sult  of  an  explosion.  The  thrust  faulting  hypothesis  has  the  added  advantage  that  the  induced  fault¬ 
ing  is  much  smaller  in  fault  dimensions,  as  measured  in  terms  of  seismic  moment,  for  thrust 
rather  than  strike-slip  to  explain  the  same  Rayleigh-wave  observations.  The  largest  moment  that 
we  obtained  for  the  Eastern  Kazakh  events  we  have  studied,  sufficient  to  reverse  phase  at  all  azi¬ 
muths,  is  quite  small  (2  x  10*33  dyn-cm).  Assuming  a  shear  modulus  of  3  X  10*1  dyn/cm3,  this 
is,  for  example,  equivalent  to  slip  of  10  cm  over  a  fault  of  dimensions  2x3  km. 
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This  choice  between  strike -slip  and  thrust  as  the  mode  of  Induced  faulting  could  be  resolved 
by  the  large  difference  in  Love -wave  excitation  between  the  two  faulting  types.  Unfortunately, 
this  information  cannot  be  deduced  from  the  data  given  by  Toksoz  and  Kehrer,  and  we  do  not  have 
available  high-quality  digital  data  of  SRO  type  with  good  azimuthal  coverage  for  NTS.  A  lengthy 
analysis  of  LRSM  data  for  earlier  events  would  possibly  resolve  the  question  of  whether  induced 
faulting  at  NTS  is  of  strike -slip  type,  as  assumed  by  Toksoz  and  Kehrer,  or  of  thrust  type,  sup¬ 
porting  the  generality  of  the  model  of  Masse. 

A  somewhat  cruder  method  of  confirming  one  or  the  other  mode  of  induced  faulting  may  be 
provided  by  the  variation  of  Mg  with  yield.  Figure  III-4  shows  the  radiation  pattern  produced 
by  the  addition  of  faulting  of  strike-slip,  thrust,  and  normal  type  to  the  explosion  monopole,  the 
amount  of  faulting  being  described  by  the  F-factor  —  here  defined  as  the  ratio  of  the  moment  of 
the  double  couple  representing  the  faulting  to  that  of  the  monopole.  From  this,  we  may  see  that 
the  radiation  pattern  produced  by  the  addition  of  induced  strike-slip  faulting  with  T  =  3.0  can  be 
reproduced  in  shape  and  phase  by  induced  thrust  faulting  of  much  smaller  moment  (ir  ~  0.7).  In 
the  former  case  the  amplitudes  are  larger  than  that  for  monopole  only,  and  in  the  latter  they  are 
smaller.  From  the  radiation  patterns  shown  in  Fig.  III-4,  we  calculated  the  effect  upon  aver¬ 
age  M  of  induced  faulting  of  strike -slip,  thrust,  and  normal  type  as  a  function  of  F,  the  ratio 
of  the  fault  moment  to  that  of  the  explosion  monopole.  Average  Mg  effect  has  been  determined 
from  the  algorithm  of  the  ratio  of  the  azimuthally  averaged  amplitude  for  (monopole  +  fault)  to 
that  for  the  monopole  only,  and  the  results  are  shown  in  Fig.  III-5.  We  see  that  the  addition  of 
normal  faulting  sharply  increases  Mg,  the  addition  of  strike-slip  initially  has  no  net  effect  then 
increases  M  ,  and  that  the  addition  of  thrust  faulting  initially  causes  a  large  decrease  (up  to 
0.5  unit)  in  Mg  and  then  increases  it  as  F  grows  larger.  An  alarming  feature  of  induced  normal 
faulting  is  that  it  would  be  difficult  to  resolve  from  surface-wave  data,  as  it  can  never  cause 
telltale  phase  reversals  and  simply  enhances  amplitudes  in  two  quadrants.  If  normal  faulting 

can  be  induced  by  explosions  it  can  substantially  increase  M  ,  and  since  normal  faulting  regions 

s 

often  involve  high  heat  flow  and  other  indications  of  low  Q,  the  P-wave  amplitudes  may  be 
smaller  than  usual.  A  combination  of  the  two  effects  may  defeat  the  M^m^  discriminant. 

The  effects  of  induced  faulting  shown  in  Fig.  Ill  -  5  are  for  average  Mg.  For  NTS  it  is  diffi¬ 
cult,  however,  to  obtain  the  good  azimuthal  coverage  at  teleseismic  distances  to  average  the 

radiation  pattern.  The  radiation  patterns  and  station  data  shown  in  the  figures  of  Toksoz  and 
14 

Kehrer  indicate  that,  in  fact,  the  observing  stations  are  concentrated  in  the  azimuths  of  one 
of  the  smaller  lobes  of  the  radiation  pattern.  In  the  direction  of  the  largest  lobes  (NW  and  SE), 
however,  the  M  effect  is  particularly  accentuated,  as  seen  in  Fig.  III-4.  In  these  directions 
the  addition  of  strike-slip  (F  =  3.0)  enhances  Mg  by  0.7  unit,  while  the  addition  of  thrust  to  pro¬ 
duce  the  same  radiation  pattern  in  shape  and  phase  (F  ~  0.7)  reduces  Mg  by  0.2  unit.  Such  a 
difference  and,  in  particular,  the  large  enhancement  of  Mg  by  strike-slip  faulting  should  be 
observable  in  these  directions. 

We  have  accordingly  measured  Mg  at  WWSSN  station  COL  (College,  Alaska),  which  the 
figures  of  Toksoz  and  Kehrer  show  to  be  in  the  direction  of  the  larger  lobe  for  events  with  large 
deduced  F-values,  for  13  events  of  known  yield  in  the  range  56  to  1200  kT.  These  events  in¬ 
clude  7  for  which  Toksoz  and  Kehrer  deduced  F-factors  assuming  strike-slip  faulting,  2  for 
which  they  had  high  F-factors  (PUedrlver,  56  kT,  F  =  3.2;  and  Greeley,  825  kT,  F  =  1.6)  and 
phase  reversals  at  some  azimuths.  The  other  6  events  are  included  to  provide  more  calibration 
of  Mg  as  a  function  of  yield;  some  of  these  may  possibly  show,  on  closer  examination,  induced 
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YIELDS,  SURFACE-WAVE  MGNITUDE,  AND  F-FACTORS 
FOR  UNDERGROUND  NUCLEAR  EXPLOSIONS  IN  NEVADA 


Event 

Explosion 

!  .. 

Y  (kT) 

M 

s 

F 

1 

Duryea 

65 

3.64 

0.75 

2 

Chartreuse 

70 

3.75 

0.90 

3 

Piledriver 

56 

3.56 

3.2 

4 

Half  Beak 

300 

4.64 

0.67 

5 

Greeley 

825 

4.83 

1.6 

6 

Scotch 

150 

4.32 

- 

7 

Knickerbocker 

71 

3.90 

- 

8 

Boxcar 

1200 

5.03 

- 

9 

Benham 

1100 

5.08 

- 

10 

Flask 

107 

3.51 

- 

11 

Carpetbag 

220 

4.02 

- 

12 

Miniata 

80 

3.60 

- 

13 

Starwort 

85 

3.51 

— 

Yields  Y  (kT)  From  Springer  and  Kinnaman,1^  in  which  dates, 
locations,  and  shot  medium  are  given  for  these  events.  Ms  as 
measured  at  WWSSN  station  COL  (College,  Alaska).  F-fdctor 
as  determined  by  Toksoz  and  Kehrer,^  assuming  induced  strike- 
slip  faulting. 


faulting  effects.  Table  III -3  lists  the  events  used,  their  announced  yield/  and,  where  deter¬ 
mined,  the  F-factors  of  Toksoz  and  Kehrer.  Their  dates,  shot  times,  depths,  and  locations 

15 

may  be  found  in  the  tables  of  Springer  and  Kinnaman. 

Figure  III -6  shows  the  M  values  measured  at  COL  as  a  function  of  yield.  Also  shown  is 

s 

the  relationship 


Mg  =  log  Y  +  2.0 


(yield  Y  in  kT) 


determined  by  Marshall  et  al.  for  explosions  in  consolidated  rock.  They  found  that,  for  events 

in  unconsolidated  rock,  M  values  were  up  to  1  magnitude  unit  smaller.  We  see  that  the  above 

formula  fits  all  the  events  used  here  to  better  than  0.5  magnitude  unit:  in  the  cases  where  our 

events  are  the  same  as  those  of  Marshall  et  ah,  the  measurements  made  at  COL  agree  very  well 

with  their  average  Mg  values.  Both  our  results  and  those  given  in  Fig.  1  of  Marshall  et  al.  show 

no  anomalously  large  M  for  Piledriver  and  Greeley,  as  would  be  expected  by  the  addition  of 

s 

large  amounts  of  strike-slip  faulting. 

Toksoz  and  Kehrer  cite  as  confirmation  of  the  induced  strike-slip  faulting  that  they  propose 
the  agreement  of  many  of  the  fault  strikes  they  obtained  with  the  strike  of  the  nearby  Yucca  fault. 
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Unfortunately  for  their  case,  moat  observed  displacements  on  this  fault  have  been  dip-slip  in 
character  and  they  were  forced  to  propose  that  observations  of  en-echelon  fractures  across  this 
fault  after  one  underground  explosion  (Corduroy)  indicated  strike-slip  movement  in  the  under¬ 
lying  basement.  Their  Figs.  8  and  9  show,  however,  very  large  numbers  of  dip-slip  faults  of 
precisely  the  orientation  required  of  thrust  faulting  to  enhance  amplitudes  in  the  NW-SE  quadrants 
and  to  reduce  them  (and  if  strong  enough,  to  reverse  the  phase)  in  the  two  other  quadrants.  We 
concede  that,  in  view  of  the  well-known  extensional  character  of  this  region,  the  purely  tectonic 
motion  along  these  faults  is  more  likely  to  have  normal,  rather  than  thrust  faulting,  type;  but  it 
is  not  inconceivable  that  this  original  direction  of  fault  movement  can  be  reversed  by  the  com- 
pressional  motion  involved  during  an  explosion. 

We  thus  believe  that  the  observations  of  azimuthal  amplitude  variations  and  phase  reversals 
for  some  events  at  NTS  are  more  consistent  with  induced  thrust,  rathpr  than  strike -slip,  faulting. 
This  implies  that  the  explosion  source  model  of  Masse,  involving  explosion  monopole  and  induced 
thrust  faulting,  with  the  possible  addition  of  spall,  is  substantially  correct  at  least  in  its  predic¬ 
tion  of  the  nature  of  surface-wave  radiation  from  underground  nuclear  explosions.  Further  con¬ 
firmation  should  be  obtained  by  an  analysis  of  Love -wave  data  from  NTS  explosions. 

However,  we  do  not  agree  with  Masse  that  the  thrust  faulting  component  does  not  imply 
triggering  of  an  earthquake.  The  consistency  of  the  fault  strikes  obtained  by  tnis  study  for 
Eastern  Kazakh  events,  and  by  Toksoz  and  Kehrer  for  NTS,  indicates  that  the  direction  of  a 
local  or  regional  tectonic  stress  field  governs  that  of  the  induced  faulting.  Symmetric  thrust 
faulting  on  a  cone  above  the  explosion  hypocenter,  equivalent  to  a  compensated  linear  vector 
dipole  source,  produces  Rayleigh  waves  reversed  in  initial  phase  compared  with  the  explosion 
itself.  If  sufficiently  large,  this  conical  faulting  may  reverse  the  phase  at  all  azimuths:  it  can¬ 
not,  however,  produce  the  enhanced  radiation  in  some  directions,  and  phase  reversals  at  some 


azimuths  only,  that  are  seen  for  both  NTS  and  Eastern  Kazakh  explosions.  Phase  reversals  at 
all  azimuths  have  not,  to  our  knowledge,  been  observed  for  NTS  events  -  this  may  be  due  to  the 
low  competency  of  the  surficial  layers  of  volcanically  derived  sediments,  that  are  pervasive  in 
the  Basin  and  Range  province  of  the  Western  U.S.,  being  unable  to  support  or  transmit  large 


local  or  regional  stresses. 


R.G.  North 


T.  J.  Fitch 


C.  CRUST  AND  UPPER-MANTLE  STRUCTURE  OF  THE  MIDDLE  EAST 

AND  SOUTH  CENTRAL  ASIA 

In  this  study,  we  invert  the  available  Rayleigh-wave  phase-  and  group-velocity  data  to  de¬ 
termine  crust  and  upper-mantle  structure  in  the  Middle  East  and  South  Central  Asia. 

Figure  III-7  shows  some  of  the  paths  for  which  velocities  are  available.  The  data  come 
from  our  measurements  (Tubman  and  Toksoz17)  as  well  as  those  of  other  investigators. 

Figures  III-8  to  III-10  show  phase  and  group  velocities  for  Iran,  the  Arabian  Peninsula,  and 
Tibet,  respectively.  The  path  in  Iran  is  between  stations  SHI  and  MSH.  Arabian  Peninsula  ve¬ 
locities  are  for  the  path  Red  Sea  -  Sin.  Velocities  for  Tibet  are  from  Bird,18  Patton  19  and  Chen 
and  Molnar.^9 

Phase  and  group  velocities  were  inverted  simultaneously  in  order  to  determine  crust  and 
upper-mantle  structure.  A  maximum -likelihood  method  weighted  in  both  data  and  model  space 
was  used.  Inversions  were  done  for  shear  velocities  only.  Compressional  velocities,  densities, 
and  layer  thicknesses  were  held  constant  during  the  iterations.  They  were  adjusted  separately 
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when  starting  further  iterations  with  a  new  starting  model.  Iterations  were  continued  until  the 

theoretical  dispersion  curves  matched  the  observed  velocities  to  within  a  specified  rms  ■.* ror. 

The  structures  for  Arabia,  Iran,  and  Tibet  are  shown  in  Fig.  Ill- 11.  The  model  for  Tibet 

21 

is  based  on  a  refraction  study  which  displays  a  low -velocity  zone  in  the  middle  crust.  The 
surface-wave  data  do  not  have  sufficient  resolution  to  indicate  the  presence  of  this  low-velocity 
layer. 

Iran  and  Tibet  exhibit  thickened  crusts  and  low  velocities.  Both  regions  have  been  deformed 
extensively  by  collisional  processes.  The  Arabian  Peninsula  is  on  the  other  side  of  the  colli- 
sional  zone  from  Iran  and  remains  relatively  unaffected.  Crustal  thicknesses  determined  are 
38,  47,  and  70  km  for  Arabia,  Iran,  and  Tibet,  respectively. 

Figures  III-8  to  III- 1 0  show  the  fit  of  the  theoretical  dispersion  curves  predicted  by  these 
models  to  the  data.  As  can  be  seen  in  the  figures,  theoretical  and  observed  velocities  agree 

qulte  weU*  K.  M.  Tubman 
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D.  NUMERICALLY  STABLE  CALCULATION  OF  MODE  EIGENFUNCTIONS 

The  solution  to  the  numerical  problems  associated  with  evaluating  the  dispersion  function 

22 

for  P-SV  in  a  layered  medium  has  been  elaborated  by  Abo-Zena  and  computational  improve- 
23  24 

ments  given  by  Menke  and  Harvey.  These  methods  avoid  the  evaluation  of  propagator  prod¬ 
ucts  that  can  have  exponentially  large  terms  that  algebraically  cancel  but  cannot  be  success¬ 
fully  machine -canceled  without  loss  of  precision.  Synthesis  of  seismograms  for  a  buried  source, 
however,  entails  not  only  evaluation  of  the  dispersion  function  but  also  evaluation  of  mode  eigen¬ 
functions.  This  process  again  encounters  the  numerical  problem  associated  with  evaluating 

24 

products  of  propagators.  Harvey  introduces  pseudolayers  and  the  constraint  imposed  by  an 
exponentially  decaying  boundary  condition  to  avoid  this  numerical  problem.  This  method  still 
relies  on  machine  cancellation  of  exponentially  positive  terms,  but  the  pseudolayers  are  kept 
sufficiently  thin  that  these  terms  are  small. 

The  exponentially  decaying  boundary  condition,  however,  can  be  used  to  construct  an  algo¬ 
rithm  for  eigenfunction  calculation  that  completely  avoids  the  use  of  pseudolayers.  This  algo- 

25 

rithm  preserves  the  advantage  of  the  asymptotic  fundamental  matrix  method  given  by  Cormier 
in  minimizing  the  number  of  vertically  inhomogeneous  layers  needed  to  describe  an  Earth  model. 
The  outline  of  the  computational  technique  that  follows  can  be  applied  to  a  homogeneously  layered 
model  as  well  as  an  inhomogeneously  layered  model. 

The  stress -displacement  eigenfunctions  at  radius  r._1  are  related  to  those  at  radius  r.  by 
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where  P  is  the  propagator  matrix  from  r._j  to  r..  Following  Harvey,  by  using  Abo-Zena's 
Y  matrix  elements  and  the  Sommerfeld  radiation  condition,  the  two  traction  eigenfunctions  can 
be  expressed  in  terms  of  the  two  displacement  eigenfunctions: 
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e3  =  Y14//Y34ei  ~  Y24  /Y34e2 

e4  =  Y13  ^Y34el  +Y23/Y34e2 


(HI-2) 


where  the  Y  matrix  elements  and  eigenfunctions  are  evaluated  at  the  same  depth.  Substituting 

i  - 1  1  - 1 

Eq.  (III-2)  into  Eq.  (Ill- 1 )  shows  that  the  eigenfunctions  e^  ,  e2  at  radius  r^  are  related  to 
e|,e2  at  r^  by  a  2  x  2  propagator  matrix  E: 


(III-3) 


where  the  elements  of  E  are  given  by 


E1l=<PllV34-P13Y14tP14Y13>',Y34 
E!2  *  <P32Y34  -  PnY34  *  PI4Y23>/Y34 
E21-(P21Y34-P23Y14*P24Y13>''Y34 
E22  =  *P22Y34  +  P23Y24  +  P24^Y34 


(III-4) 


Y  elements  being  evaluated  at  radius  r^,  and  P  being  the  propagator  from  r.  _  ^  to  r..  F6r  an 
Earth  model  having  n  —  i  discontinuities  between  radius  r^  and  the  free  surface,  Eq.  (III-3)  can 
be  generalized  by 


ElE1+1 .  .  .  En 


(III-5) 


where  the  index  e  is  the  ellipticity  at  the  free  surface,  and  E  is  the  eigenfunction  propagator 
from  the  first  model  boundary  below  the  free  surface  to  the  free  surface.  At  any  depth  the  trac¬ 
tion  eigenfunctions  e^,  can  be  found  by  Eq.  (Ill -3).  By  substituting  specific  expressions  for  the 
Y  elements  and  propagator  elements  in  Eq.  (III-4)  and  exploiting  all  the  algebraic  cancellations 
that  occur,  it  can  be  shown  that  the  numerator  of  each  element  of  an  E1  matrix  consists  of  four 
terms  of  different  exponential  order.  The  order  of  each  of  these  terms  is  always  less  than  or 
equal  to  the  exponential  order  of  the  denominator  YJ4.  When  each  numerator  term  is  normalized 
by  the  same  factor  as  Y^,  numerical  stability  is  insured  without  need  of  inserting  pseudolayer 

boundaries.  ,  r  ,,  „  . 

V.  F.  Cormier 

E.  THE  WAVEGUIDE  OF  Lg 

The  regional  seismic  phase  termed  "Lg"  refers  to  surface  shear  waves  having  a  dominant 
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period  0.5  to  6  s  and  an  average  group  velocity  of  3.5  km/s  (Press  and  Ewing  ).  Lg  can  be  ob¬ 
served  on  all  three  components  of  ground  motion,  although  it  is  commonly  largest  on  the  hori- 
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zontal  component  (e.g.,  see  Street  et  aL  ).  Because  Lg  is  often  the  largest  arrival  at  regional 

distances,  seismic  discriminants  based  on  Lg  have  particular  value  for  smaller  events. 

28 

Gupta  et  aL  tested  an  earthquake -explosion  discriminant  formed  from  the  ratio  of  the  maximum 
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P-wave  amplitude  to  the  maximum  Lg  amplitude  (P  /Lg  ).  Chen  and  Pomeroy  tested 

ili  HX  III  cuv 

a  discriminant  given  by  the  ratio  of  the  energies  of  Lg  In  two  group-velocity  windows  (E^/Ej, 

where  E^  Is  the  energy  in  the  window  4.0  to  3.4  km/s  and  Ej  is  the  energy  in  the  window  3.4  to 

2.8  km/s).  Given  that  the  properties  of  Lg  can  be  successfully  predicted  from  a  superposition 

of  higher-mode  Rayleigh  and  Love  waves  in  the  appropriate  group-velocity  window  (Oliver  and 
30  31  32 

Ewing  ;  Knopoff  et  aL  ;  Panza  and  Calcagntle  ),  a  study  of  the  waveguide  of  higher  surface - 
wave  modes  combined  with  an  examination  of  their  group-velocity  dispersion  can  be  used  to 
optimize  discriminants  based  on  Lg  amplitudes.  One  result  of  such  a  study  would  be  a  predic¬ 
tion  of  the  group-velocity  and  frequency  windows  of  Lg  that  are  best  excited  by  a  deep  source  vs 
those  windows  best  excited  by  a  near-surface  source. 

As  a  first  step  in  such  a  study,  the  relative  amplitudes  and  eigenfunction  depth  behavior  of 
the  higher  Rayleigh  modes  were  computed  for  a  crustal  model.  The  computational  techniques 

used  are  described  in  the  preceding  section  of  this  report.  Figure  111-12  shows  the  group  ve- 
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locities  of  Rayleigh  modes  of  the  Southern  California  model  of  Kanamori  and  Hadley  in  the 
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window  appropriate  for  the  Lg  phase.  (This  crustal  model  is  same  as  that  used  by  Harvey  in 

a  test  of  seismogram  synthesis.  The  model  calculated  here,  however,  includes  the  effects  of 

Earth  sphericity  by  the  use  of  asymptotic  wave  functions  as  discussed  in  Cormier.25)  In 

Fig.  111-12  one  can  see  narrow  bands  of  group  velocity  in  which  the  stationary  phases  of  groups 

of  higher  modes  cluster.  This  clustering  of  stationary  phases  may  account  for  the  observations 

of  multiple  Lg  phases  described  by  Bath.  For  the  discussion  that  follows,  the  modes  in  the 

frequency  window  1.3  to  2  Hz  will  be  referred  to  as  Lg3  in  the  group-velocity  window  3.35  to 

3.50  km/-c,  Lg,  in  the  window  3.24  to  3.32  km/s,  and  Lg,  in  the  window  2.96  to  3.06  km/s. 

“  35  J 

Panza  et  aL  define  the  normalized  energy  density  of  a  mode  in  terms  of  the  density  func¬ 
tion  p(z),  the  vertical  and  horizontal  displacement  functions  u(z),  w(z),  and  their  values  at  the 
free  surface: 

E  =  p(z)  {[u(z)/u(o)]2  +  (w(z)/w(o)]2}  .  (III-6) 

In  terms  of  ellipticity  and  displacement  eigenfunctions,  the  energy  density  is  given  by 

E  =p(z)  [e2  +  (e2/e)2]  .  (III-7) 

The  frequencies  at  which  the  mode  energy  densities  of  Lgj,  Lg2,  and  Lg3  were  calculated  were 

chosen  to  be  at  or  near  points  of  stationary  phase  in  the  mode  group -velocity  curves.  In  each 

energy-depth  plot  (Figs.  Ill- 1 3  to  III-15),  the  S-velocity  profile  is  shown  for  comparison. 

Figure  III- 13  shows  the  energy  of  a  mode  member  of  Lgj  to  be  maximum  in  the  depth  range 

5  to  30  km.  Mode  members  of  Lg2  (Fig.  III-14)  and  Lg3  (Fig.  Ill- 1 5),  in  contrast,  have  large 

energy  throughout  the  depth  range  0  to  35  km  and  maximum  energy  in  the  uppermost  crustal 

layer.  From  these  results  it  can  be  concluded  that  the  crustal  layers  form  the  waveguide  of 

Lg  -  the  free  surface  forming  the  upper  boundary  of  the  waveguide,  and  the  moho  forming  the 

lower  boundary.  The  thinness  of  oceanic  crust  implies  that  if  an  oceanic  phase  analogous  to  Lg 

does  exist,  it  will  have  a  characteristic  frequency  much  higher  than  that  of  continental  Lg.  For 

continental  models  having  a  sedimentary  layer,  the  earliest  portion  of  the  Lg  wavetrain  in  the 

highest  group-velocity  window  has  most  of  its  energy  concentrated  in  the  crustal  layer  lying 

beneath  the  sedimentary  layer.  The  waveguide  of  this  portion  of  the  Lg  wavetrain  is  consistent 
26 

with  Press  and  Ewing's  idea  that  Lg  propagates  in  the  granitic  layer  of  the  crust. 
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The  waveguide  of  Lg^  la  also  consistent  with  the  success  of  the  seismic  discriminant  con¬ 
structed  by  Chen  and  Pomeroy  using  the  ratio  of  Lg  energies  in  different  group-velocity  windows. 
That  study  found  E^/Ej  1°  be  smallest  for  the  Salmon  explosion  compared  with  Ef)/Ej  for  earth¬ 
quakes  in  eastern  North  America.  From  Fig.  111-13,  It  is  apparent  that  an  explosion  at  or  near 
the  free  surface  of  a  sedimentary  layer  will  be  less  efficient  In  exciting  the  modes  that  Eh  mea¬ 
sures  compared  with  an  earthquake  beneath  the  sedimentary  layer.  If  this  interpretation  is 
correct,  the  E^/E^  discriminant  would  not  perform  as  well  for  Lg  propagating  in  a  crust  lacking 
a  sedimentary  layer.  These  results  suggest  the  feasibility  of  using  the  energy  depth  behavior  of 
Lg  modes  to  optimise  seismic  discriminants.  Once  the  average  crustal  structure  that  charac¬ 
terizes  an  Lg  path  is  known,  optimal  frequency  and  group-velocity  filters  can  be  constructed  to 
discriminate  source  depth.  v>  F<  Cormier 

F.  INVERSION  FOR  MULTIPLE  SOURCES 

4 

This  section  describes  generalization  of  the  method  of  centroid  location  *  to  the  problem 
of  searching  for  the  source  parameters  (moment  tensor,  hypocentral  location)  of  either  inde¬ 
pendent  events  closely  spaced  in  time  or  events  of  finite  spatial  and  temporal  dimensions.  In 
particular,  we  shall  demonstrate  application  of  this  method  in  a  numerical  experiment  in  which 
the  mechanism  of  a  strike-slip  earthquake  and  explosion  of  commensurate  moments  are  fully 
resolved.  The  following  describes  the  principle  of  the  approach. 

Assume  that  there  are  N  point  sources  at  starting  locations  (r  ,  6  ,  <p  )  and  origin  times 

t  .  In  the  starting  iteration,  an  attempt  is  made  to  determine  their  moment  tensors  f '  for 
n  a  in 

the  k  record  of  our  data  set,  the  equation  of  condition  is: 


N  G 

uk(r,t)  =2  E  ff 
n=l  i«l 


(III -8) 


After  the  appropriate  system  of  normal  equations  is  solved  for  the  initial  estimates  of  f.'^\ 
the  elements  of  the  moment  tensor  and  location  parameters  are  refined  in  successive  iterations: 

uk(£*t)  ~uk°^r»t)  *  E  an°)(t)  6rn  +  bn°>(t)  6Gn +  cn°>(t)  6^n 


*  dr«>  “n  + 


(III-9) 


The  appropriate  definitions  and  derivatives  are  given  for  a  single  source  in  Ref.  36. 

In  Table  III-4  we  give  a  synthetic  example  of  application  of  this  approach  to  the  resolution 
of  two  sources  approximately  150  km  apart  and  activated  within  10  s  of  each  other.  The  first 
source  was  a  strike-slip  earthquake  (f2  =  1,  f^  =  -1;  see  Ref.  1  or  36  for  explanation  of  notation); 
the  second  was  an  explosion  (f^  -  =  f^  =  1).  The  starting  locations  are  quite  remote  from  the 

true  ones,  and  the  initial  iteration  yields  moment  tensors  that  have  little  in  common  with  the 
true  source  mechanism  of  the  individual  events.  However,  after  several  iterations,  both  sources 
are  fully  resolved  and  their  locations,  origin  times,  as  well  as  the  moment  tensors  agree  very 
well  with  the  exact  values.  Although  this  has  been  a  synthetic  example,  the  result  is  of  some  in¬ 
terest,  as  we  are  not  aware  of  any  other  method  that  would  allow  such  resolution  with  a  compara¬ 
bly  sparse  receiver  network. 
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TABLE  1 1 1-4 

RESOLUTION  OF  A  DOUBLE  SOURCE* 


Iteration  0 


Iteration  I 


Iteration  2 


iteration  3 


Iteration  4 


S.1  S.2  S.l  S.2 


Latitude  35.50  37.50  36.02  36.50 

Longitude  77.50  79.50  78.42  78.81 

6tn  -3.00  6.00  0.28  11.22  I  1.31 


36.02  37.00  36.00 

78.06  79.01  78.01 

-0.26  10.43  0.02 


-0.561  0,963  0.371  0.436  0.034  0.947  -  0.009  0.983  -  0  <V,2  0.997 

0.517  1.016  1.296  0.508  1.027  0.975  0.992  0.990  0.999  1.002 

-1.542  1.256  -  0.271  0.140  -  0.943  0.928  -  0.997  0.989  -1.001  0.998 

0.738  -  0.911  -1.007  0,755  -  0.066  -  0.032  -  0.012  0.008  -  0.002  0.001 

0.559  -  0.565  -  0.627  0.493  -  0.043  -  0.005  -  0.003  -  0.002  0.000  -  0,002 

0.191  -0.106  -  0.169  0,123  -  0.031  0.001  -0.003  -  0.002  0.001  0.000 


*An  experiment  with  synthetic  data.  Seismograms  were  computed  for  a  network  of  13  SRO/ASRO  sta¬ 
tions  assuming  the  following  source  parameters:  (1)  Strike-slip  earthquake  (f2  =  1,  fg  *  - 1,  f]  =  f4  = 
fj  =  f^  =  0)  at  location  36°N,  78®E,  and  origin  time  0  S;  and  (2)  an  explosion  (f ^  =  f 2  ~  ^3  =  ^  ^  ^4  = 
f^  =  f^  =  0)  at  location  37°N,  79°E,  and  origin  time  10  s.  The  sum  of  these  two  sets  of  seismograms 
represented  the  data  set  in  inversion.  The  starting  coordinates  are  given  in  ihe  column  for  Iteration  0 
(only  moment  tensor  is  obtained  in  the  zeroth  Iteration).  S.l  and  S.2  represent  sources  1  and  2, 
respectively. 


Another  potential  application  of  this  approach  is  in  studies  of  sources  of  finite  size.  Imagine 
that  an  earthquake  can  be  approximated  by  a  propagating  line  source,  with  the  fault  length  L  and 
duration  2T.  The  centroid  method"^  for  a  single  source  would  yield  a  point  at  the  center  of  the 
fault  and  shift  in  the  origin  time  equal  to  T.  As  discussed  in  the  next  section  (Sec.  G),  this  in¬ 
formation  is  sufficient  to  resolve  the  fault  parameters,  provided  that  neither  the  location  of  the 
point  of  origin  nor  the  centroid  location  are  biased  by  the  effect  of  lateral  heterogeneities.  Our 
current  experience  indicates  that  such  bias  is  probable,  and  thus  the  estimates  of  the  fault  length 
velocity  and  direction  of  rupture  propagation  could  be  grossly  erroneous.  If  the  propagating 
fault  were  to  be  represented  by  two  point  sources,  their  locations,  ideally,  should  fall  on  the 
fault  at  points  distant  by  1/4  and  3/4L  from  the  point  of  origin,  and  at  times  1/2  and  3/2 T  after 
the  origin  time.  However,  even  if  due  to  the  large-scale  heterogeneities  the  absolute  locations 
are  biased,  in  this  case  the  estimate  of  the  source  dimensions  should  be  much  more  accurate,  as 
the  ray-paths  from  both  point  sources  to  the  receiver  are  very  similar. 

Even  though  our  current  experience  with  tile  analysis  of  actual  data  is  limited,  the  initial 
results  are  encouraging,  or  at  least  intriguing.  The  main  difficulty  appears  to  be  associated  with 
assuring  the  stability  of  the  solution  by  imposing  physically  plausible  smoothing  conditions. 

A.  M.  Dziewonski 


G.  ROUTINE  DETERMINATION  OF  SOURCE  PARAMETERS  FROM  THE  ANALYSIS 
OF  WAVEFORMS  RECORDED  BY  THE  GLOBAL  DIGITAL  NETWORK 

In  1980,  Chou  and  Dziewonski*  proposed  a  method  of  derivation  of  the  elements  of  the  seis¬ 
mic  moment  tensor  through  linear  inversion  of  the  body-wave  sections  of  records  of  the  SRO/ 
ASRO  digital  seismograph  network.  The  fundamental  mode  data  are  excluded  from  the  analysis, 
because  of  their  sensitivity  to  the  lateral  variations  in  the  Earth's  structure.  Derivation  by 
Woodhouse  of  the  differential  kernels  for  perturbation  of  synthetic  seismograms  with  respect 

to  changes  in  the  geographical  coordinates,  depth,  and  the  origin  time  of  an  event  allowed  joint 
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inversion  for  the  source  mechanism  and  location  of'the  centroid  of  the  stress  glut.  The  de¬ 
tailed  description  of  the  procedure  with  numerous  examples  of  application  to  studies  of  global 

36 

and  regional  seismicity  has  been  given  by  Dziewonski,  Chou,  and  Woodhouse.  Since  that  report 

has  been  completed,  we  have  analyzed  approximately  100  events  ranging  in  seismic  moment  from 

2x10  to  2  x  10  dyn-cm.  The  ability  to  cover  this  wide  a  spread  of  earthquake  sizes  using 

exactly  the  same  procedure  is  the  result  of  the  wide  dynamic  range  of  the  SRO/ASRO  network. 

It  is  possible  to  think  of  extending  the  current  lower  limit  of  the  seismic  moment  to 
23 

2  x  10  dj.:-cm  (Mg  ~  4.5)  by  calibrating  the  source -receiver  paths  for  dispersion  and  attenua¬ 
tion  of  the  fundamental-mode  surface  waves.  This  requires  occurrence  of  earthquakes  with 
Mg  ~  6  in  a  given  region.  The  suggested  procedure  is  as  follows.  Once  the  source  mechanism 
is  determined  for  a  large  event,  using  the  approach  of  Ref.  36,  the  phase  velocity  and  apparent 
attenuation  can  be  determined  for  the  Rayleigh-wave  train  contained  in  the  same  set  of  records. 
These  calibration  curves  can  be  used  to  determine  moments  and  source  mechanism  of  consider¬ 
ably  smaller  events  in  the  vicinity  of  the  master  earthquake.  Substantial  effort  will  be  necessary 
to  establish  the  range  of  distances  from  the  master  event  over  which  the  calibration  curves  are 
valid.  Most  likely,  there  will  be  strong  dependence  on  the  tectonic  nature  of  the  source-receiver 
path;  one  would  anticipate  that  paths  along  continental  margins,  mid-ocean  ridges,  or  trenches 
could  be  rather  unstable. 

The  next  subject  to  be  discussed  in  this  report  has  to  do  with  rather  large -magnitude  (6.5  to 

8)  events.  Because  of  the  response  of  the  SRO/ASRO  network,  the  energy  of  the  signals  used  in 

our  analysis  is  strongly  concentrated  near  a  60-s  period  (the  spectrum  is  tapered  and  truncated 

at  a  45-s  period).  For  the  events  of  magnitude  above  7,  or  events  with  anomalous  time  function 

("slow"  earthquakes),  it  is  possible  that  our  results  could  seriously  underestimate  the  seismic 

moment.  For  this  reason,  we  have  expanded  our  proc  >dure  to  incorporate  information  contained 
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in  records  of  the  ultra-long  period  IDA  (International  Deployment  of  Accelerometers  )  network. 
The  dynamic  range  of  these  instruments  is  rather  small  (2048  digital  counts),  and  following  a 
major  earthquake  the  initial  portion  of  the  record  is  distorted.  However,  these  records  contain 
excellent -quality  information  on  excitation  of  waves  with  periods  greater  than  100  s  which,  de¬ 
pending  on  the  Q  of  the  given  mode,  may  be  detected  from  hours  to  months  (only  qSq,  really) 
after  an  event.  In  reaching  the  decision  how  long  a  section  of  record  to  use,  one  should  recog¬ 
nize  the  fact  that  our  knowledge  of  the  elastic  and  anelastic  structure  of  the  Earth  is  imperfect. 

If,  "or  a  norma)  mode,  the  difference  between  the  actual  and  computed  frequency  and  attenua¬ 
tion  factors  of  a  particular  normal  mode  is  <5wk  and  So^,  the  difference  between  predicted  and  ob¬ 
served  displacements  is: 

“akt 

<5uk  ~  t(6a>k  sinw^t  +  6c»k  cos  w^t)  e 
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Thua,  for  small  perturbations,  the  difference  increases  linearly  with  time.  The  conclusion  is 
that  one  should  use  time  series  beginning  as  close  to  the  origin  time  as  possible,  and  minimize 
their  length,  providing  that  all  pertinent  information  about  the  source  properties  has  been  in¬ 
cluded.  In  case  of  IDA  data,  this  minimum  length  is  approximately  3  h:  the  time  required  by 
the  slowest  mantle  waves  to  travel  once  around  the  world.  With  an  unimportant  exception  of  the 
modes  that  correspond  to  nearly  vertical  ScS  reflections,  all  other  modes  propagate  with  higher 
group  velocities,  and  we  are  assured  that  their  contributions  from  both  minor  and  major  arc 
propagation  will  be  included  in  a  3-h-long  segment  of  a  record.  It  is  worthwhile  to  note  that 
quite  different  strategy  would  be  advantageous  if  the  purpose  of  data  analysis  were  to  measure 
functionals  of  the  Earth's  structure  such  as  periods  of  normal  modes  and  their  attenuation. 

Another  important  decision  is  related  to  the  problem  of  scaling  of  the  information  contained 
in  the  data  from  two  sets  of  instruments  of  very  different  response.  Decision  has  been  made  to 
weigh  the  records  from  the  IDA  network  in  such  a  way  that,  on  average,  contribution  to  the  sys¬ 
tem  of  normal  equations  of  the  data  from  an  IDA  record  would  be  the  same  as  that  of  data  from 
a  record  from  an  ASRO  or  SRO  station. 

Figure  III  - 16  shows  distribution  of  SRO/ASRO  and  IDA  stations  that  were  used  in  studying 
the  Irplnia  earthquake  of  23  November  1980.  This  magnitude  7  earthquake  in  southern  Italy 
caused  death  of  at  least  3,000  persons  and  destroyed  living  quarters  of  300,000  people.  The  map 
is  an  equidistant  azimuthal  projection  centered  on  the  epicenter.  The  short -broken  lines  connect 
to  the  epicenter  the  SRO  (squares),  ASRO  (up-triangles),  and  DWWSSN  stations  (down-triangle). 
The  long-broken  lines  are  full  great-circles  passing  through  an  appropriate  IDA  station,  as,  in 
this  case,  information  on  the  waves  excited  in  both  minor  and  major  arc  directions  is  included. 

Figures  III-17(a)  through  (c)  show  examples  of  comparison  between  the  observed  records 
(upper  trace)  and  synthetics  computed  for  our  solution  for  the  source  parameters.  The  param¬ 
eter  AMAX  indicates  the  SRO  digital  count  for  the  maximum  amplitude  for  a  given  pair  of  sta¬ 
tions;  the  response  of  DWWSSN  and  ASRO  stations  has  been  appropriately  normalized. 

The  observed  and  computed  traces  for  the  IDA  records  are  shown  in  Fig.  111-18.  The 
parameter  DELAY  represents  the  time  interval  (in  hours)  between  the  origin  time  and  beginning 
of  a  particular  record  shown  in  the  figure.  The  agreement  is  very  good.  The  dispersion  of  sur¬ 
face  waves  in  the  selected  records  does  not  exhibit  major  differences  with  respect  to  the  refer¬ 
ence  model  (1066B  (see  Ref.  39)],  but  cases  of  major  shifts  have  been  encountered  on  other 
occasions. 

The  very  good  match  of  both  SRO/ASRO  as  well  as  the  long -period  IDA  data  has  been 
achieved,  in  part,  by  introducing  a  finite  duration  of  the  source:  half-time  duration  of  15  s. 

The  apparent  shifts  in  the  origin  time  of  17  s  are  very  large  for  an  event  of  this  magnitude.  This 
indicates  the  possibility  that  this  might  have  been  a  multiple  event. 

Table  III- 5  summarizes  our  results  of  inversion  for  the  Irpinia  earthquake.  Also  included 
in  the  table  are  results  for  the  Yugoslavian  earthquake  of  15  April  1979.  This  earthquake  had 
a  seismic  moment  very  similar  to  the  Irpinia  earthquake,  but  its  mechanism  was  entirely  dif¬ 
ferent.  While  the  Irpinia  event  was  a  normal  fault,  the  Yugoslavian  earthquake  was  a  thrust 
event.  The  strikes  of  both  earthquakes  are  nearly  parallel.  One  might  speculate  that  there  is 
a  link  between  these  events  on  the  opposite  sides  of  the  Adriatic  micro-plate.  Figures  III-19(a) 
through  (c)  and  III -20  show  comparison  between  the  data  and  synthetics  for  the  event  of 


TABLE  1 1 1-5 

SOURCE  PARAMETERS  OF  TWO  EARTHQUAKES* 


23  November  1980 

15  April  1979 

Origin  Time 

18:34:53.3 

6:19:44.1 

Latitude  (°N) 

40.89 

42.10 

Longitude  (°E) 

15.33 

19.21 

StQ  (*) 

17.1 

16.6 

Latitude  (°N) 

40.73 

42.12 

Longitude  (°E) 

14.97 

18.68 

Depth  (km) 

25 

(10) 

Half-Duration  (s) 

15 

15 

fl 

-3.05  ±0.05 

2.75  ±  0.05 

f2 

1.08  ±0.07 

-1.32  ±0.06 

f3 

1.97  ±  0.04 

-1.43  ±0.05 

f4 

0.96  ±0.16 

1.24  ±0.14 

f5 

-0.94±  0.16 

-  1 .16  ±  0.17 

,f6 

-0.48  ±  0.06 

0.77  ±  0.06 

Moment 

2.47 

3.28 

Plunge/AZM 

13°/62° 

73°/41° 

Moment 

0.92 

-0.61 

Plunge/AZM 

6°/331° 

1  °/136° 

Moment 

-3.39 

-2.68 

[  Plunge/AZM 

76°/218° 

17  0/226° 

*  Relocation  parameters  and  the  elements  of  the  moment  tensor  for  the  Irpinia 
earthquake  of  23  November  1980  and  Yugoslavian  earthquake  of  15  April  1979. 
Note  that  the  two  earthquakes  are  of  similar  size  but  the  elements  of  their 
moment  tensors  are,  essentially,  opposite  in  sign.  The  elements  of  the  moment 
tensors  and  the  principal  values  are  multiplied  by  a  scale  factor  of  lO^dyn-cm. 
The  focal  depth  of  the  Yugoslavian  earthquake  was  held  fixed  at  10  km. 
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It  is  Interesting  to  note  that  the  relative  contribution  of  overtones  is  different  for  the  Italian 
(depth  25  km)  and  Yugoslavian  (depth  10  km)  earthquakes.  Even  though  the  difference  in  depth 
is  relatively  small,  the  actual  recordings,  as  well  as  synthetics,  show  substantial  increase  in 
the  relative  amplitude  of  overtones  for  the  Italian  earthquake. 

One  interesting  conclusion  of  the  analysis  of  these  two  large,  damaging  earthquakes  is  that 
highly  reliable  information  on  the  source  properties  can  be  obtained  from  a  relatively  sparse 
network.  If,  in  the  future,  data  from  globally  distributed  digital  stations  were  transmitted  in 
real  time,  it  is  possible,  I  think,  to  obtain  information  on  the  size  and  mechanism  of  important 
earthquakes  within  a  few  hours  after  an  event.  In  cases  of  tsunamogenic  earthquakes,  this  could 
be  used  for  issuing  appropriate  warnings.  For  earthquakes  occurring  near  densely  populated 
areas,  estimates  of  the  earthquake  size,  direction,  and  length  of  the  fault  could  be  used  for  early 


estimates  of  potential  damage  and  needed  aid. 


A.  M.  Dziewonski 


H.  EXACT  NUMERICAL  SOLUTION  FOR  A  PROPAGATING  LINE  SOURCE 
AND  THE  METHOD  OF  CENTROID  LOCATION 


The  most  commonly  used  mathematical  model  of  an  earthquake  of  finite  spatial  and  temporal 
dimensions  is  a  propagating  line  source.  This  approximation  is  justified  for  earthquakes  with 
the  fault  lengths  significantly  exceeding  their  widths.  In  addition  to  its  practical  applications, 
the  model  is  attractive  from  a  theoretical  point  of  view,  because  for  a  radially  (vertically)  sym¬ 
metric  Earth  model  variation  of  displacements  in  the  horizontal  plane  is  described  by  analytical 
functions.  Thus,  it  is  possible  to  derive  relatively  simple  expressions  for  excitation  of  surface 
waves  or  free  oscillations  generated  by  such  a  source.  On  the  other  hand,  numerical  integration 
is  necessary  to  consider  the  effect  of  source  propagation  in  the  vertical  plane  and  it  is  difficult, 
if  not  impossible,  to  examine  the  functional  dependence  between  the  source  parameters  and  radia¬ 
tion  of  seismic  waves. 

Generation  of  surface  waves  by  a  propagating  line  source  was  first  considered  by  Hen- 
40 

Menahem,  who  showed  that  the  spectrum  of  a  wavetrain  due  to  a  point  source  is  modified  by 
the  directivity  factor: 


sinX  _-iX 
X  e 


(III-10) 


where 


v  Lw  , C  . 

X  =  2C  “  ~  cosy) 


(III-ll) 


and  L  is  the  fault  length,  is  the  frequency,  C  is  phase  velocity,  v  is  the  velocity  of  rupture 
propagation,  and  y  is  the  angle  between  the  direction  of  rupture  propagation  and  direction  to 
the  receiver. 

This  expression  has  been  utilized  in  numerous  studies,  most  notably  by  Kanamori,  who  first 
applied  it  to  the  analysis  of  the  excitation  of  mantle  waves  generated  by  the  Kurile  Islands  earth¬ 
quake  of  1963  (Ref.  41)  and  the  Alaskan  earthquake  of  1964  (Ref.  42),  and  subsequently  to  records 
of  many  major  events  of  this  century. 

As  the  expression  (III- 10)  applies  to  contribution  of  a  wavegroup  for  a  particular  dispersion 
branch  (most  frequently  fundamental  modes)  and  a  specific  orbit  (Rj,  for  example),  it  is  difficult 
to  use  it  in  the  analysis  of  waveforms  corresponding  to  superposition  of  numerous,  strongly  dis¬ 
persed  overtone  branches.  This  used  to  impose  limitations  on  the  type  of  data  that  could  be  ap¬ 
plied  in  studies  of  the  effect  of  finiteness  of  the  source. 
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In  1977,  Dziewonski  ai'd  Romanowlcz  derived  an  exact  solution  to  the  problem  of  excitation 
of  normal  modes  by  horizontal  line  source  representing  a  fragment  of  a  great  circle.  If  the  co¬ 
ordinate  system  is  rotated  such  that  the  fault  lies  in  the  equatorial  plane  with  the  direction  of 

rupture  propagation  coincident  with  the  increasing  longitude,  then  the  expressions  of  Gilbert  and 
39 

Dziewonski  for  the  strain  at  the  source  can  be  integrated  analytically  along  the  fault  length. 
There  is  a  simple  relationship  between  excitation  of  a  normal  mode  due  to  a  point  source  and  the 
line  source.  The  displacement  speptrum  of  the  kth  mode  of  angular  order  number  l  excited  by 
the  point  source  is: 

uk(r,  w)  =  hk(w  )  Yi  3™^™  (111*12) 

m--i 

■  a  orkt 

where  hk(co)  is  the  spectrum  of  h^ft)  =  1  —  e  cosw^t,  and 


G 


i=l 


(III-13) 


with  f.  being  the  six  independent  elements  of  the  moment  tensor,  and  the  strains  e  defined  by 
1  39 

Gilbert  and  Dziewonski  in  Eqs.  (2.1.19)  for  spheroidal  modes  and  (2.1.20)  for  toroidal  modes; 
the  vector  s  is  defined  in  Eq.  (2.1.1)  of  that  paper.  Dziewonski  and  Romanowicz  have  shown 
that  the  effect  of  the  line  source  is  to  modify  the  excitation  parameters  ip™  by  a  frequency- 
dependent  factor: 


where 


A  m 

$  Jtw) 


sin  X.. 


m 


-iX 


m  ,  m 
*k 


1  lur o  \ 
X  =  7  *  — -  +  mj 
m  2  \  v  / 


(III— 44) 


(III- 15) 


and  4*  is  the  angular  extent  of  the  fault,  rQ  is  the  source  radius,  :■  id  v  is  the  velocity  of  rup¬ 
ture  propagation.  They  have  also  pointed  out  the  similarity  bet*-  jen  the  directivity  terms  in 
Eqs.  (III-10)  and  (III-12). 

While  Eqs.  (Ill - 12 )  and  (III-14)  provide  the  means  for  the  exa -t  evaluation  of  the  spectra  of 
normal  modes  excited  by  a  line  source,  the  frequency-domain  representation  would  be  very 
inefficient  in  the  calculation  of  synthetic  seismograms  in  applications  considered  by  Dziewonski, 
Chou,  and  Woodhouse,3^  who  add  contributions  from  roughly  5000  normal  modes  to  obtain  a 
wavetrain  composed  of  long-period  body  waves. 

It  may  be  noted  that  the  spectrum 


sinX  -iX 
_ m  p  m 

m 


is  the  Fourier  transform  of  a  function  f(t): 


for  t  <  0 


f(t)  =  Jr  exPl~imvt/r0J 


for  0<  t-$  2T 


for  t  >  2T 


where  T  =  (1/2)  4>rQ/v  is  the  half -duration  of  the  source. 

Thus,  to  account  for  the  finite  size  of  the  source,  one  needs  to  evaluate  the  convolution 
integral: 


■»»  ■  A  y 

O 


^  pmin(t,  2T)  -imvr/ro 


-ok(t-r) 


<(t -r)] 


(III— 16 ) 


This  integral  may  be  readily  evaluated  analytically.  For  reasons  of  numerical  efficiency, 
it  is  expedient  to  consider  excitation  at  times  greater  than  2T  (that  is,  after  the  fault  has 
stopped),  disregard  contributions  of  permanent  displacements,  approximate  through  lineariza¬ 
tion  the  effect  of  attenuation  during  the  process  of  fault  propagation,  and  to  change  the  limits  of 
integration  to  (— T,  T).  The  result  is: 


Sm(t-T>=e 


°k<t  f)  f(A™  +  B™)  coswk(t  -  T)  +  i(A™  -  Bkm)  sincok(t  -  T)] 


where 


sin  (uk  +  0m)  T 


2(0).  +/3  T 
4  k 


3in(o)k-/im)T 

2(a)k-/3m)T 


with  3m  3  mv/V 


The  displacement  vector  associated  with  the  k  normal  mode  is  thus: 
-ak(t-T) 

uk(r,t-T)  =  e  K  [Ck  coswk(t-T) +  Sj.  sinwk(t -T)] 


(III- 17) 


where 


£k  =  Re  2  <5m§km<Akm  +  Bkm> 


(III -18) 


where  is  the  Newman's  delta  (6  =1  for  m  =  0;  6  =2  for  m  >  0)  and  ib,  are  the  excitation 

m  m  m  k 

coefficients  of  Eq.  (Ill- 1 3 )  for  a  point  source  in  the  "equatorial"  coordinate  system.  The  time 
shift  of  T  cam  be  applied  to  the  final  synthetic  seismogram.  The  advantage  of  the  procedure 
implied  by  Eqs.  (III-17)  and  (III- 18 )  over  that  associated  with  Eq.  (Ill- 14)  is  that,  in  the  former 
case,  it  is  necessary  to  apply  summation  to  <1  +1)  scalar  values  while,  in  the  latter  case,  it 
would  be  necessary  to  add  (f  +  1)  spectra,  each  represented  by  several  hundred  values. 
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The  execution  time  needed  for  computation  of  a  set  of  synthetic  seismograms  for  a  propa¬ 
gating  line  source  assuming  typical  coverage  of  an  event  by  12  SRO/ASRO  three -component 
stations  and  10  IDA  vertical  component  gravimeters  is  about  25  min.  on  an  Interdata  8/32  com¬ 
puter,  some  four  to  five  times  greater  than  for  the  point  source  —  a  substantial  but  not  formi¬ 
dable  difference. 

Equation  (111-16)  may  be  modified  to  evaluate  excitation  for  source  time  functions  other 
than  a  boxcar.  In  general: 

pmin(t,2T)  -Imvr/r  t  -a.(t-T) 

gm(t)  =  J  a(r)  e  |l-e  coswk(t-r) 

where  the  source  function  a(r)  is  normalized  such  that: 


|  dT  (111-19) 


a(r)  dT  =  1 


There  is  a  wide  class  of  functions  appropriate  for  representation  of  the  source  (i.e.,  triangle, 
trapeze)  for  which  integral  (III- 19)  can  be  evaluated  analytically. 

It  should  also  be  pointed  out  that,  for  a  given  set  of  parameters  related  to  the  rupture  propa¬ 
gation,  excitation  of  normal  modes  remains  linear  with  respect  to  the  zero-frequency  components 
of  the  moment  tensor.  Thus,  the  equations  presented  here  can  be  readily  used,  after  obvious 
modifications,  in  formulation  of  an  inverse  problem  for  the  moment  tensor  of  a  propagating  line 
source. 

The  solid  lines  in  Fig.  III-21  represent  a  set  of  synthetic  seismograms  for  a  lone  source  of 
a  length  of  1°  of  arc  (111.2  km)  propagating  due  east  with  a  velocity  of  3.3  km/s.  The  source 
mechanism  is  a  45°  normal  fault  at  a  depth  of  15  km  with  the  strike  along  the  east -west  axis. 

All  receivers  spaced  with  a  45"  step  in  azimuth  are  located  at  a  distance  of  90"  and  have  response 
of  the  SRO  instruments.  The  dotted  lines  are  the  difference  between  the  result  for  propagating 
line  source  and  the  point  source  at  the  centroid  location,  corrected  for  the  finite  duration.  The 
effect  of  the  direction  of  source  propagation  on  amplitudes  is  clearly  discernible.  However,  the 
zero-crossings  of  the  difference  traces  coincide  with  those  of  the  signal.  Thus,  a  point  source 
at  the  centroid  location  accounts  very  well  for  the  phase  shifts  associated  with  the  propagating 
line  source.  This  is  explained  by  the  following  consideration. 

The  change  in  the  travel  time  of  a  wave  with  a  phase  velocity  C  from  the  point  of  origin  to 
a  receiver  at  a  distance  x  is  the  sum  of  the  perturbation  in  the  origin  time:  6t  =  L/v;  and  the 
change  in  the  travel  time  due  to  the  perturbation  in  the  epicentral  distance.  To  the  first  order: 

6x  =  —(1/2)  L  cosy,  where  y  is  the  angle  between  direction  of  rupture  propagation  and  direction 
to  the  receiver.  Thus,  the  total  change  in  the  travel  time  is 

6t  =  fe  (v  ~  coay) 


and  the  phase  shift  at  a  frequency  u : 

6<p  =  w6t  =  -  cosy) 


an  expression  identical  with  Eq.  (III-ll). 

This  indicates  that  the  method  of  centroid  location^  can  be  used,  in  principle,  to  resolve 
the  highly  nonlinear  problem  of  the  propagating  line  source,  that  is  to  determine  the  fault  length. 
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velocity,  and  direction  of  rupture  propagation.  This  conclusion  has  been  verified  by  experiments 
on  synthetic  data  using  realistic  distribution  of  receivers.  However,  in  processing  actual  seis¬ 
mograms  one  should  be  aware  that  the  eplcentral  location  as  well  as  the  parameters  of  the  cen¬ 
troid  may  be  biased  by  the  effect  of  lateral  heterogeneities. 

Perhaps  the  effects  of  mislocation  and  the  source  size  could  be  resolved  by  studying  many 
earthquakes  of  a  wide  range  of  magnitudes  from  the  same  source  region,  thuB  developing  an 
equivalent  master  event  technique.  Another  possible  approach  is  outlined  in  Sec.  G  above. 

A.  M.  Dziewonski 


I.  THE  EXCITATION  OF  SEISMIC  WAVES  BY  A  SOURCE 
SPANNING  A  SURFACE  OF  DISCONTINUITY 

The  displacement  seismograms  due  to  any  seismic  source  may  be  written  as  a  function  of 
spatial  coordinates  x  and  time  t  (see  Refs.  37.  39,  and  44) 


s(x,t)  =  Yj  2k(x)  hk(t)  *  *k(t)  (III— 2 0 ) 

k 

where  the  summation  is  over  all  normal  modes,  and  where  *  denotes  convolution.  The  vector 
fields  sk<x)  are  the  normal  mode  eigenfunctions,  and  the  time  functions  (t)  are  given  by 

-cv  t 

lyt)  =  1-e  K  coswkt  (HI —21) 

ij. 

where  wk  and  o:k  are  the  eigenfrequency  and  attenuation  constant  of  the  kl  mode.  The  excitation 
coefficients  ^k(t)  in  Eq.  (III.-20)  are 


(III —22 ) 


k  th  * 

where  e  is  the  elastic  strain  in  the  k  mode  (overbar  indicates  complex  conjugation),  and  T  is 

the  glut  rate  distribution  of  the  source.  The  spatial  integration  in  Eq.  (III-22)  is  over  the  source 

region  Vg. 

This  representation  of  the  seismograms,  in  terms  of  a  superposition  of  normal  modes,  is, 

of  course,  under  suitable  approximations,  equivalent  to  representations  more  often  used  in 

45 

body-wave  seismology,  as  has  recently  been  explicitly  shown  by  Chapman  and  Woodhouse.  The 
results  to  be  presented  here  are  generally  applicable  and  are  not  restricted  to  the  normal-mode 
representation. 

If  the  source  region  V  is  sufficiently  small  compared  with  the  wavelengths  of  modes  included 
s  k 

in  Eq.  (III-20),  it  will  be  an  adequate  approximation  to  assume  e  (x)  to  be  constant  in  Vg,  and 
Eq.  ( III  -  22 )  then  gives  the  point  source  approximation 

^k(t)  =  ?k(x0)  :  M(t)  (III-23) 


where  Xq  is  a  point  in  Vg,  and  M(t)  is  the  moment  rate  tensor  of  the  source 
M(t)  =  (  f(x',t)  d3x'  . 


(III —24 ) 
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It  is  the  purpose  of  this  contribution  to  point  out  that  this  argument  fails  if  the  source  re¬ 
gion  V  spans  a  surface  of  discontinuity,  such  as  the  moho,  and  to  investigate  the  effect  of  such 
a  circumstance  on  the  excitation  of  seismic  waves  and  upon  the  interpretation  of  moment  tensors 
obtained  from  inversion.  The  results  presented  below  show  that  the  seismic  radiation  from  a 
source  spanning  a  discontinuity  is  equivalent,  in  the  point  source  approximation,  to  that  for  a 
particular  source  placed  on  one  side  of  the  discontinuity,  or  to  that  of  a  different  source  placed 
on  the  other  side,  and  that  significantly  different  moment  tensors  can  be  obtained  from  inversion, 
depending  upon  which  side  of  the  discontinuity  the  source  is  assumed  to  lie.  Neither  of  these 
equivalent  sources  corresponds  to  the  moment  rate  tensor  defined  in  Eq.  (111-24),  and  it  is  not 
possible,  in  general,  to  recover  the  true  moment  rate  tensor  unless  further  assumptions  about 
the  source  are  made.  A  set  of  assumptions,  which  seem  reasonable  for  certain  physical  situa¬ 
tions  and  which  make  the  inverse  problem  determinate,  is  suggested  below.  An  important  con¬ 
sequence  of  this  indeterminacy  is  that  constraints,  such  as  that  the  moment  tensor  be  trace  free 
or  that  it  represent  a  double-couple,  cannot  be  applied  unless  one  makes  such  further  assump¬ 
tions,  since  the  constraints  will  have  different  effects,  depending  upon  which  side  of  the  discon¬ 
tinuity  one  assumes  the  source  to  lie. 

Consider  a  small  source  region  Vs  which  is  divided  into  two  subregions  V+,  V_  by  the  sur¬ 
face  2,  across  which  the  elastic  moduli  are  discontinuous.  We  shall  define  a  local  Cartesian 

coordinate  system  in  which  the  x,  direction  is  normal  to  2  and  denote  the  six  independent  com¬ 
ic  3q  1 

ponents  of  e  and  M(t)  by 

el  =  ell  '  e2  =  e22  *  e3  =  e33  ’  e4  =  ei2  *  e5  =  e13  *  e6  =  e23 

fl  =  M11  *  f2  =  M22  *  f3  =  M33  *  f4  =  M12  *  f5  =  W13  ’  f6  =  M23 
In  place  of  Eq.  (Ill -2  3)  we  now  have 

d'k(t)  =  e+  •  f+  +  e"  •  f"  (111-25) 

where  e*  are  the  strains  evaluated  on  each  side  of  2 ,  and  f*  are  source  vectors  corresponding 
to  the  moment  rate  tensors 


r(x'.t)  d3x' 


(III-26) 


Because  of  the  required  continuity  of  the  modal  displacements  and  tractions  across  2,  e±  must 
satisfy  certain  continuity  conditions;  the  tangential  strains  e2<e3*e6  must  be  continuous,  together 
with  the  normal  tractions  (X  +  2^^  +  X(e2  +  e3),  Spe^,  2p,e5,  where  X,p,  are  the  Lame  param¬ 
eters.  It  follows  that  e*  in  Eq.  (III-25)  satisfy 


A+e+  =  A'e" 

where  A*  are  the  matrices  with  nonvanishing  elements 


A11  =  x±  +  z^±  * 

A*  -  A  *  - 
A22  “  A33  ' 

A4*4  =  A55  =  Z»*±  ■ 

A  *  -  a  * 
A12  "  A13 
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Equation  (III-25)  may  now  be  written  in  such  a  way  that  it  depends  only  upon  e+  or  only 
upon  e“,  i.e.. 


(t)  3  e+  •  F+  a  e*  •  F* 


(111-2?) 


where 


F+  3  f+  +  (A  1A+)T  f" 
F"  3  f"  +  <A+'1A')Tf+ 


(111-28) 


where  superscript  T  indicates  matrix  transposition.  It  is  clear  from  Eq.  (111-27)  that  the  true 
source  is  equivalent  to  the  source  F+  placed  on  the  plus  side  of  £  or  to  the  source  F*  placed 
on  the  minus  side  of  £.  and  F"  represent  the  two  different  sources  which  would  be  obtained 
from  inversion  of  seismic  data  under  the  assumption  that  the  source  lay  on  the  plus  or  the  minus 
side  of  the  discontinuity  surface  £.  Neither  of  these  sources  represents  the  true  moment  rate 
tensor  components 

f  3  f+  +  f  " 

corresponding  to  the  moment  tensor  defined  in  Eq.  (111—24).  From  Eq.  (III-27),  we  find  that  the 

■f 

equivalent  sources  F  ,  F  are  related  to  each  other  by 


Specifically, 


Fi  =AF1 


.A^'-V  +  A-^V 


(III -29) 


f2+  =  f;-bf; 


F3  =F3  ~BF1 


F4  =CF4 


F5  =  CF5 


F6+  =  Fb 


(IH-30) 


where 


A  =  (A+  +  2pJ/(A_  +  2  n_) 
B  =  (A_  -A+)/(A.  +  2n_) 

c  ~  P+/V_  • 


(III-31 ) 


For  a  major  discontinuity,  suc.i  as  the  moho,  or  an  interface  between  slow  sedimentary  strata 
and  basement  rocks,  the  sources  F+,  F-  can  differ  quite  appreciably,  both  in  magnitude  and  in 
geometry.  Using,  for  instance,  the  parameters  from  the  model  1066B2  appropriate  for  the  moho, 
and  taking  T  to  indicate  above  the  moho,  we  find 
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A  a  0.545  ,  B  a  0.02? 


C  a  0.434 


If  two  solutions  are  obtained  for  the  moment  tensor  of  a  seismic  source,  one  assuming  It  to  lie 
just  above  the  moho  and  the  other  assuming  it  to  lie  just  below,  they  are  only  directly  comparable 
if  transformed  by  the  above  relations.  Furthermore,  if  F1-,  say,  is  constrained  to  represent  a 
source  without  volume  change  (F^+  +  F^  +  F3  =  0),  or  to  represent  a  double-couple,  then  F  will 
not  have  the  corresponding  property.  If  the  source  is  explosive,  but  is  near  a  major  discontinu¬ 
ity,  it  can  appear  to  have  a  significant  devlatoric  part. 

The  true  moment  tensor  components  f  =  f+  +  f '  cannot  be  determined  from  the  seismic  ra¬ 
diation  without  making  further  assumptions  about  the  sour  ce.  One  possible  assumption  would  be 
that  f+,  f"  are  proportional,  as  would  be  the  case,  for  example,  if  a  single  fault  cut  across  the 
discontinuity.  Let  us  assume,  then,  that 


f  =  <1  —  #3)  f 


(111-32) 


where  /)  is  a  parameter  representing  the  fraction  of  the  total  moment  lying  on  the  minus  side  of 
the  discontinuity.  In  this  case,  we  have  from  Eq.  (III-28) 


F+  =  [/?  +  /J(a"1a+)t]  L 

•F'  =  [/3  +  0(A+'V)T]f 


(III -33) 


where  (3  *  1  —  /3. 

Let  us  suppose  that.  F+  is  the  source  obtained  from  inversion  under  the  assumption  that  the 
source  lies  on  the  plus  side  of  the  discontinuity.  Then,  Eq.  (III-33)  gives 

ft  ■  (?  +  A (If1  F  +  I 


f2  =  F2  +  /3B</3  +  A/3)"1  Fj 
f3  =  F3+  +  /3B(0  +  A /3)“  F* 


f4  =  (/3C  +  0)"1  F4+ 


(III-34) 


f.  =  (/3C  +  (3)"1  F$+ 


f6  "  F6 


where  A.  B,  and  C  are  those  given  in  Eq.  (III-31).  If  it  is  required  to  constrain  f  to  represent 
a  source  without  volume  change  (f3  +  f2  +  f3  =0),  then  the  inversion  for  F+  should  be  carried  out 
under  the  corresponding  constraint 


(1  +  2/3)  F1  +  (IS  +  A/3)  <F2+  +  F3+)  =  0 


(111-35) 


Similarly,  if  f  is  to  represent  a  double-couple  the  corresponding  nonlinear  constraint  must  be 
translated  into  a  constraint  on  F+. 

With  the  above  prescription,  inversion  for  the  true  moment  tensor  f  may  be  carried  out  by 
making  only  minor  modifications  to  existing  inversion  algorithms.  The  result,  however,  will 
depend  upon  the  specified  parameter  /3.  If  a  least-squares  procedure  is  used,  these  modifica¬ 
tions  can  be  performed  during  the  final  stage  of  inversion  in  which  only  the  inner  product  matrix 


and  the  corresponding  vector  are  required;  thus,  many  values  of  /3  may  be  tried  at  little  compu¬ 


tational  cost. 


J.  H.  Woodhouse 


J.  A  LINEAR  APPROXIMATION  TO  THE  EFFECT  OF  FINITE  RUPTURE  LENGTH 
AND  RUPTURE  VELOCITY  ON  THE  RADIATION  FROM  SEISMIC  SOURCES 

In  our  ongoing  effort  to  retrieve  from  long -period  waveform  data  as  much  Information  as 
possible  about  seismic  sources.  It  Is  desirable  to  Incorporate  the  effect  of  a  finite  rupture  length 
and  rupture  velocity.  It  should  be  possible  to  obtain  direct  estimates  of  these  quantities  by  In¬ 
cluding  the  extra  parameters  In  our  inversion  procedure.  In  order  to  begin  the  linearisation  in¬ 
version  for  the  rupture  velocity,  it  is  necessary  first  to  have  a  linear  approximation  to  its  effect 
on  the  predicted  waveforms;  it  is  our  purpose  here  to  derive  such  an  approximation. 

The  extra  solution  for  a  rupture  of  infinitesimal  vertical  and  lateral,  but  of  arbitrary,  length, 
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propagating  along  a  great  circle,  has  been  given  by  Dziewonski  and  Romanowicz  and  was  dis¬ 
cussed  in  Sec.  II  above.  This  exact  solution,  however,  does  not  provide  a  convenient  starting 
point  to  obtain  the  linearized  result  sought  here. 

We  shall  consider  a  "line  source,"  viewed  in  a  Cartesian  coordinate  frame  in  which  the  rup- 
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ture  propagates  along  the  x-axis.  The  source  will  be  represented  by  the  glut  rate  distribution 

r<x,t)  =  2^  6(x  -  vt)  <5(y)  6(z)  H(t  +  T)  H(T  -t)  M  (III-36) 

where  M  is  the  total  moment  tensor,  v  is  the  rupture  velocity,  assumed  constant,  and  2T  is  the 
source  duration.  In  Eq.  (Ill- 36 ),  6(x)  and  H(x)  are  the  Dirac  6-function  and  the  Heaviside  step 
function,  respectively.  The  excitations  of  free  oscillations  for  this  source  are  obtainable  directly 
from  Eq.  (Ill -22 )  above.  On  performing  the  spatial  inteh.ations  we  find,  simply, 

57s  M  :  ek(x)  |t|  <  T 

1  x=vt,y=0,  z=0 

♦k«>  = 

0  1 1 1  >  T  (III-37) 

We  shall  now  make  the  approximation  that  the  modal  strains  e  (x)  are  sufficiently  well  approxi¬ 
mated  along  the  rupture  by  the  first  two  terms  in  their  Taylor  expansions  about  the  midpoint. 

We  then  have 

“  FT  M  :  [§k(9)  +  *  el  §k(~)  0]  1*1  <  T  •  (111-38) 

Convolving  with  the  modal  time  functions  h^t),  as  in  Eq.  (III-36),  we  find 


/  sinw.T  -o.t  \  v 

\(t)  *  ^k(t>  *  ( 1  -  „  t  e  cos  v )  m  =  r<s 


_  v  / _ 

u  \  w 


sin u.  T 
k 


sinwkt  M  :  ~  (x) 


where  terms  proportional  to  the  small  quantity  cv^  have  been  neglected. 


(111-39) 


(111-40) 


/  Blnw.T  -a.t  v  k 

lik(t)  *  =  (* - uTT- e  cos  w^t  j  M;e 

k 

/  ainw.  T  .  -«kt  v.  . 

-(  ^r--coeukT)  e  sUi^kt—M:?! 

where  v.  are  the  three  components  of  the  rupture  propagation  velocity  vector  (summation  over 

1  k  k 

1  assumed)  and  where  the  modal  strains  e  and  strain  gradients  e  ,  are  those  evaluated  at  the 

"  137  '•* 

midpoint  of  the  rupture,  which  is  also  the  source  centroid.  *  Similarly,  time  in  Eq.  (111-39) 
is  measured  relative  to  the  centroid  time,  namely  a  time  T  after  the  onset  of  rupture. 

A  successful  inversion  procedure  for  the  estimation  of  the  centroid  location  and  time  from 
waveform  data  has  been  described  in  an  earlier  issue  of  SATS1  and  the  formulas  for  the  strain 
gradients  e1^  in  Eq.  (111-40)  have  also  been  given  previously.1,36  By  virtue  of  Eq.  (111-20), 

Eq.  (111-39 )  provides  a  direct  linear  relationship  between  the  predicted  waveforms  and  the  com¬ 
ponents  of  the  rupture  velocity  vector,  and  will  form  the  basis  for  a  linearized  inversion  for 
these  quantities  together  with  an  iterative  inversion  for  the  half-duration  T  of  the  source.  This 
technique  will  be  applied  first  to  synthetic  data  calculated  from  the  exact  solution  for  a  line 
source  as  discussed  in  Sec.  H  above.  It  will  then  be  applied  to  a  combined  data  set  consisting 
of  long-period  body  waves  (~16  MHz)  from  SRO  and  ASRO  records  and  long-period  surface 
waves  (~6  MHz)  in  IDA  records,  for  a  number  of  earthquakes. 

J.  H.  Woodhouse 

K.  THE  CALCULATION  OF  THE  COMPLETE  NORMAL- MODE  SPECTRUM 
OF  THE  EARTH  USING  FINITE  DIFFERENCE  METHODS 

The  synthesis  of  complete  long-period  seismograms  by  the  superposition  of  normal  modes 

forms  the  basis  of  many  seismological  studies,  and,  in  particular,  forms  the  basis  of  some  of 
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the  studies  of  the  seismic  source  reported  here  and  in  an  earlier  SATS.  Clearly,  it  is  desir¬ 
able  to  extend  such  studies  to  higher  frequencies.  The  major  advantage  of  this  method  of  syn¬ 
thesis  is  that  it  is  completely  automatic  and  omits  no  part  of  the  theoretically  predictable  seis¬ 
mic  signal.  To  use  the  method,  however,  one  must  have  previously  calculated  a  complete 
catalog  of  normal-mode  eigenfrequencies  and  eigenfunctions,  up  to  some  preassigned  limit  in 

frequency.  At  present,  only  two  complete  catalogs  are  available,  comprising  the  modal  data 
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up  to  22  MHz  for  the  models  1066A  and  1066B,  calculated  by  Buland  using  a  Rayleigh-Ritz 
procedure. 

It  is,  of  course,  possible  to  calculate  normal -mode  eigenfrequencies  and  eigenfunctions  by 
solving  numerically  the  radial  differential  equations,  but  this  method  has  been  incapable  of  pro¬ 
ducing  a  complete  catalog  of  modes  in  a  reasonable  computation  time  since  no  systematic  way 
of  searching  for  roots  of  the  secular  equation  has  been  devised.  Simply  to  scan  all  frequencies 
requires  such  a  small  frequency  step  as  to  make  the  method  quite  impracticable.  The  major  ad¬ 
vantage  of  the  Rayleigh-Ritz  technique  is  that  it  converts  the  problem  into  a  matrix  eigenvalue 
problem,  for  which  a  Sturm  count  can  be  calculated.  Thus,  given  any  two  frequencies  it  is  pos¬ 
sible  to  say,  without  an  exhaustive  search,  how  many  roots  of  the  secular  equation  lie  between 
those  frequencies.  This  allows  all  the  roots  to  be  bracketed  very  efficiently,  and  convergence 
on  the  bracketed  roots  is  then  relatively  straightforward. 
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Here,  we  report  that  a  similar  scheme  is  now  available  for  the  differential-equation  pro¬ 
cedure  (if  the  gravitational  perturbation  due  to  the  oscillations  is  neglected)  since,  as  with  the 
Rayleigh-Ritz  procedure,  we  are  now  able  to  calculate  how  many  roots  of  the  secular  equation 
lie  between  two  given  frequencies  simply  by  integrating  twice  the  system  of  ordinary  differential 
equations.  It  is  well  known  that  for  second-order  Sturm-Liouville  boundary-value  problems, 
such  a  count  of  eigenvalues  may  be  established  in  terms  of  the  number  of  aero  crossings  of  some 
combination  of  the  eigenfunction  and  its  derivative.  This  method  is  applicable  to  toroidal  and 
radial  modes,  but  in  these  cases  the  modes  are  sufficiently  well  and  predictably  separated  in 
frequency  that  it  need  not  be  used.  Until  now  there  has  been  no  corresponding  result  for  spher¬ 
oidal  modes,  which  are  often  extremely  closely  and  unevenly  spaced  in  frequency. 

In  addition  to  the  difficulty  of  finding  all  roots  of  the  secular  equation,  a  further  difficulty 
in  using  the  differential-equation  approach  is  an  essential  numerical  instability  in  the  calculation 
of  the  secular  determinant.  This  difficulty  has  been  solved  by  several  authors  in  different  con¬ 
texts22, the  general  solution  of  the  problem  is  to  perform  calculations  in  terms  of  the 
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"minor  vectors"  of  pairs  (or  triplets)  of  solutions.  A  technique  for  calculating  eigenfunctions 

SI 

entirely  in  terms  of  such  solutions  of  the  minor  equations  has  recently  been  given  by  Woodhouse. 

The  technique  to  be  described  here  has,  as  yet,  been  developed  only  for  the  case  in  which 
the  gravitational  perturbation  due  to  the  oscillations  is  neglected.  This  approximation  is  prob¬ 
ably  entirely  adequate  for  frequencies  exceeding,  say,  10  MHz.  At  lower  frequencies,  correc¬ 
tions  may  be  applied  for  the  effects  of  self-gravitation,  or  the  modal  eigenfunctions  may  be  re- 
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diagonalized  to  account  for  the  gravitational  perturbation.  In  practice,  the  catalogs  of  Buland 

were  also  calculated  in  this  way  for  reasons  of  computational  economy.  The  new  technique  has 

been  used  to  calculate  all  spheroidal  modes  up  to  a  frequency  of  35  MHz  for  the  Preliminary 
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Reference  Earth  Model  of  Dziewonski  and  Anderson,  including  the  effects  of  transverse  isot¬ 
ropy  in  the  upper  mantle  and  the  frequency-dependent  corrections  to  elastic  moduli  due  to  at- 
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tenuation.  The  dispersion  diagram  for  angular  orders  up  to  100  is  shown  in  Fig.  III-22,  The 

calculation  was  performed  for  all  angular  orders  for  which  the  fundamental  mode  lies  below 

35  MHz  in  frequency,  namely  up  to  order  f  »  360. 

We  will  not  go  into  detail  here  about  how  the  algorithm  for  counting  eigenvalues  was  derived. 

In  fact,  a  rigorous  proof  of  the  result  is  stil’  lacking  though  we  can  show  that  it  has  several  of 

the  properties  required  of  such  a  count.  Here,  we  shall  simply  give  a  recipe  for  its  calculation. 

Neglecting  the  gravitational  perturbation,  the  equations  governing  spheroidal  eigenfunctions 
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may  be  written 

gj=A(r,w,f)b  .  (HI-41) 

In  solid  regions  of  the  Earth,  Eq.  (Ill -41)  is  a  fourth-order  system  of  ordinary  differential  equa¬ 
tions.  and  in  fluid  regions  (core  and  ocean)  it io  second-order  system.  In  Eq.  (Ill -41 ),  r  is 
radius,  w  is  angular  frequency,  and  t  is  angular  order;  and  the  components  of  b  are:  in  solid 
regions,  bj  =  U,  b2  =  (f(f  +  1)]1(/2  V,  b3  =  P,  b4  =  [l(i  +  1)J1//2S;  and  in  fluid  regions,  ^  =  U, 
b2  =  P  —  where  U,  V,  P,  and  S  are  functions  of  r  such  that  the  components  of  the  vector  dis¬ 
placement  field  in  the  Earth  are  given  by^ 


t  The  factors  of  r  in  these  equations  differ  from  the  conventions  of  Phinney  and  Burridge. 
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(Ill -42) 


rsr  =  U(r)Y  “(e,*) 
rsQ  =  V(r)  [l(i  +  1))‘1//2  e0 Yfm(Q,  <*>) 

rs0  =  —  V(r)  [f (f  +  l)]'1/2  coaeces ^*”(0,  4>) 
and  r,6,  <p  are  spherical  polar  coordinates.  The  components  of  ralial  traction  are 
rrrr  =  P(r)  Y™(B,<p) 

rrrQ  =  S(r)  [1(1  +  l)]"l/2 

rTrtp  -S(r)  [i(i  +  1)]*1/2  Y“(e.«) 

where  r  is  the  stress  tensor. 

At  the  center  of  the  Earth  are  two  linearly  independent  solutions  b(r)  which  satisfy  the  ap- 
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pro;-  .hie  conditions  of  regularity.  When  these  solutions  are  integrated  to  the  inner-outer  core 
boundary,  there  will  be  a  single  2 -vector  solution  in  the  core  satisfying  the  conditions  of  con¬ 
tinuity  of  U,  P  at  the  inner-core  boundary.  When  this  solution  is  integrated  to  the  core -mantle 
boundary,  there  will  be  two  linearly  independent  solutions  in  the  mantle  satisfying  the  same  con¬ 
ditions.  Similarly,  a  single  solution  in  the  ocean  will  satisfy  these  conditions  at  the  ocean  floor, 
and  the  secular  equation  is  then  obtained  by  requiring  that  the  normal  traction  P(r)  vanishes  at 
the  ocean  surface.  At  other  interfaces,  b(r)  is  continuous.  In  solid  regions,  let  us  define  the 
minor  vector  of  the  two  allowable  solutions  b^(r),  b2(r)  as 

ml  =  bllb22  ~b12b21 
m2  =bilD23  _b13b21 
=  bll  b24  _b14b21 

}  (111-44) 

m4  =  b12  b2S  _bi3b22 
m5  =  °12  b24  -b14b22 
m6  =  b13  b24  “  b14  b23 

where  the  second  subscript  is  used  to  specify  individual  elements  of  b. ,  b_. 
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It  may  be  shown  that 

m^lr)  +  m,.(r)  =  0  (III -45 ) 

and  m(r)  satisfies  a  6th-order  system  of  equations  of  the  form 
dm 

-gjr  -  A(r,w,f)m  (III -46) 
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which,  by  virtue  of  Eq.  (III-45),  immediately  reduces  *h  order.  It  may  be  shown  '  that  the 
eigenvalue  problem  may  be  stably  solved  by  integral  ly  Eq.  (III-46)  in  solid  regions  of  the 
Earth,  subject  to  the  requirement  of  continuity  at  welde.  iterfaces  and  to  suitable  boundary  con  ¬ 
ditions  at  fluid-solid  interfaces.  Or.  passing  from  a  solid  to  a  fluid,  these  are 


(III-43) 
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bl  =  m3 


b2  =m6 


(III -47 ) 


and  on  passing  from  a  fluid  to  a  solid 
ml  =  bl 


m4  =“b2 


m2  =  nij  =  irij  s  =  0 


(III-48) 


ft(r) 


Let  us  now  define  the  two  functions  of  radius 

m^(r)  in  solid  regions 

in  fluid  regions 

m3(r)-.m^(r)  in  solid  regions 
in  fluid  regions 


b2(r) 


f2(r)  = 

We  can  now  define  the  function 


bt(r) 


n(w,i)  =  £  sgn[f2(r.)  f^(r.)}  +  oc  + 


(III-49) 


(III  -  50 ) 


(III-51) 


where  the  summation  is  over  all  zeros  of  f.  (r),  and  sgn  takes  values  ±1  according  to  the  sign  of 

1  -  +' 

its  argument.  oc  is  defined  to  ne  0  or  1  dependent  on  ^  (b),  ^  (b)  being  positive  or  negative, 
where  b  is  the  radius  of  the  outer  core  and  superscripts  plus  and  minus  indicate  that  f^  or  fj1 
is  evaluated  just  outside  or  just  inside  the  core.  The  additional  quantity  <jf  in  Eq.  (111-51)  is 
simply  unity  if  i  =  1,  and  zero  otherwise,  and  is  included  simply  to  make  this  equation  concor¬ 
dant  with  the  conventional  designation  of  overtone  number. 

It  is  our  contention,  well  validated  in  numerical  tests  and  indicated  by  theoretical  considera¬ 
tions,  that  n(w,f)  is  the  spectral  function  for  the  spheroidal  mode  spectrum.  That  is  to  say, 
n(w,()  is  equal  to  the  overtone  number  of  the  mode  of  largest  frequency  less  than  w.  Conse¬ 
quently,  the  number  of  eigenfrequencies  between  two  frequencies  (u>2  >  is  given  by 

nfu^)  -  n(w^).  It  is  this  property  which  is  invaluable  in  enabling  all  eigenperiods  to  be  rapidly 

bounded,  and,  consequently,  a  complete  catalog  to  be  constructed. 

52 

In  the  calculations  for  the  PREM  model,  approximately  1.3  evaluations  of  n(w,j()  per  mode 
were  necessary  to  bracket  all  roots.  A  straightforward  procedure  of  linear  interpolation  and 
bisection  was  then  employed,  using  the  values  of  the  secular  determinant  to  converge  on  the 
roots.  This  is  the  more-time-consuming  part  of  the  process.  The  calculation  of  some 
9000  modes  took  approximately  50  h  using  an  Interdata  8/32  computer.  A  part  of  the  dispersion 
diagram  is  shown  in  Fig.  111-22.  j.  Woodhouse 
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Fig.  III-l.  Data  used  in  inversion  scheme.  Real  part  of  spectrum  integrated 
over  period  range  20  to  45  s  after  correction  back  to  source,  as  a  function  of 
station  azimuth.  Numbers  correspond  to  events  given  in  Table  III— 1. 


Fig.  III-2(a-b).  Two  comparisons 
of  short-period  P-wave  amplitudes 
with  deduced  monopole  component 
of  source  (see  text). 
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Fig.  III-3.  Observed  azimuthally  averaged 
Love -wave  spectral  amplitudes  compared 
with  that  predicted  by  deduced  thrust  fault¬ 
ing  component  of  source. 
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Fig.  Ill— 4.  Rayleigh-wave  patterns  resulting  from  addition  of  various  fault  types 
of  moment  Q  =  FE  to  a  monopole  of  strength  E,  for  different  values  of  F.  Dashed 
circles  denote  radiation  from  monopole  only.  Shaded  portions  have  phase  re¬ 
versed  compared  with  monopole. 
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Fig.  Ill- 5.  Azimuthaliy  averaged  effect  upon 
Ms  (Mg  units)  caused  by  addition  of  strike- 
slip,  normal,  and  thrust  faulting  to  a  mon- 
o'  ale  of  moment  E,  as  a  function  of  fault 
n  ament  Q  =  FE. 


Fig.  Ill -6 .  Ms  measured  at  COL  (College, 
Alaska)  as  a  function  of  yield  for  13  NTS 
explosions. 
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Fig.  1II-8.  Observed  phane  velocities  (+)  and  group  velocities  (o) 
i'or  Iran.  Lines  are  theoretical  velocities  predicted  by  derived 
model. 
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Fig.  Ill- 1 1.  Velocity  profiles  for  Iran,  Arabia,  and  Tibet. 
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Fig.  III-12.  Group-velocity  dispersion  of  Lg  modes  of  a  Southern  California 
crustal  model.  Rayleigh-mode  numbers  are  given  for  each  curve. 
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Fig.  HI  - 1 5.  Energy  density  for  Rayleigh-mode  19  at  1.6875  Hz 


Fig.  III-16.  Digital  network  coverage  of  Irpinia  event  of  23  November  1980.  Squares 
are  SRO  stations;  upward -pointing  triangles,  ASRO  stations;  downward-pointing  tri¬ 
angle,  first  operational  DWWSSN  station  (Longmire,  Washington);  and  diamonds,  IDA 
gravimeters.  Both  SRO  and  IDA  instruments  are  operational  on  Guam. 
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Fig.  lll-17(a-c).  Comparison  of  observed  (top  trace)  and  synthetic  seismograms  from 
SRO/ASRO/DWWSSN  networks  for  Irpinia  earthquake  of  23  November  1980;  see 
Table  III- 5  for  source  parameters  used  to  construct  the  synthetic  seismograms.  Pa¬ 
rameter  AMAX  is  maximum  amplitude  for  a  given  pair  of  traces  in  terms  of  digital 
counts  of  an  SRO  receiver. 
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Fig.  III-17(a-c).  Continued. 
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Fig.  Ill — 18*  Comparison  of  observed  (top  trace)  and  synthetic  seismograms  from 
IDA  network  for  Irpinia  earthquake.  Parameter  DELAY  is  time  elapsed  between 
origin  and  beginning  of  traces  shown.  Symbols  R~,  R3»  and  describe  specific 
orbits  of  mantle  Rayleigh  waves. 
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Fip.  Jir-19(a-r).  Same  as  Fig.  Ill- 17(a-e)  but  for  event  of  1  5  April  1979 
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Pig.  Ill -2 0.  Same  as  Fig.  III-18  but  for  event  of  15  April  1979.  Because  this  earthquake 
was  shallower  than  Irpinia  event,  relative  contribution  of  overtones  is  less  than  that  in 
Fig.  Ill  - 1 8 . 
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Fig.  111-21.  Solid  lines  are  synthetic  seismograms  of  long-period  body  waves 
(truncated  just  before  the  arrival  of  fundamental  modes)  for  a  propagating  line 
source  evaluated  according  to  exact  method  described  in  text.  All  modes,  except 
for  fundamentals,  in  periods  greater  than  45  s  were  used  in  calculation  (roughly 
5000  modes).  Dotted  line  represents  difference  between  exact  solution  and  a  so¬ 
lution  obtained  for  point  source  located  at  midpoint  of  fault.  Spectrum  of  point 
source  has  been  corrected  for  finite  duration.  Numbers  adjacent  to  traces  indi¬ 
cate  azimuth  of  a  receiver  in  respect  to  direction  of  fault  propagation.  See  text 
for  details. 
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G I  OSSARY 


AA 

Automatic  Association 

ABAR 

Arrival  Based  Analyst  Review 

A/D 

Analog/Digital 

ASRO 

Abbreviated  Seismic  Research  Observatory 

ATN 

Augmented  Transition  Network 

bpl 

Bits  per  Inch 

CCS 

California  Computer  Systems,  Inc. 

CPU 

Control  and  Processing  Unit 

DC 

Zero  Frequency 

DEC 

Digital  Equipment  Corporation 

DOE 

Department  of  Energy 

DP 

Detection  Processing 

DWWSSN 

Digital  World-Wide  Standard  Station  Network 

EBAR 

Event  Based  Analyst  Review 

GDSN 

Global  Digital  Seismograph  Network 

GTS 

Global  Telecommunications  System 

IDA 

International  Deployment  of  Accelerometers 

IDCE 

International  Data  Collection  Experiment 

IESD 

International  Exchange  of  Seismic  Data 

I/O 

Input -Output 

ISC 

International  Seismological  Center 

LRSM 

Long-Range  Seismic  Measurements 

LTA 

Long-Term  Average 

MARS 

Multiple -Arrival  Recognition  System 

NEIS 

National  Earthquake  Information  Service 

NSS 

National  Seismic  Station 

NTS 

Nevada  Test  Site 

PDE 

Preliminary  Determination  of  Epicenters 

RAM 

Random  Access  Memory 

rms 

Root  Mean  Square 

RST 

Remote  Seismic  Terminal 

RSTN 

Regional  Seismic  Test  Network 
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SAS 

Seismic  Analysis  Station 

SATS 

Semiannual  Technical  Summary 

SDAC 

Seismic  Data  Analysis  Center 

SDC 

Seismic  Data  Center 

SDCS 

Seismic  Data  Collection  System 

SDC -VI 

Version  1  of  Prototype  SDC 

SPZ 

Short- Period,  Vertical  Component 

SRO 

Seismic  Research  Observatory 

STA 

Short-Term  Average 

UC 

University  of  California 

uses 

U.  S.  Geological  Survey 

WMO 

WWSSN 


World  Meteorological  Organization 
World-Wide  Standard  Station  Network 
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